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Abstract 


An  analytic  and  numerical  study  of  the  problem  of  mechanical  impulse  propagation  through  a 
horizontal  alignment  of  progressively  shrinking  (tapered)  elastic  spheres  that  are  placed  between  two 
rigid  end  walls  is  investigated.  Particular  attention  is  paid  to  the  shock  absorption  and  nonlinear 
dynamical  properties  as  they  pertain  to  energy  partition.  The  studies  are  confined  to  cases  where 
initial  loading  between  the  spheres  is  zero.  The  spheres  are  assumed  to  interact  via  the  purely 
repulsive  and  strongly  nonlinear  Hertz  potential.  Two  systems  are  studied,  each  representing  a 
staggering  number  of  possible  chain  designs.  Propagation  of  energy  is  analytically  studied  in  the 
hard-sphere  approximation  and  parameter  space  diagrams  plotting  normalized  kinetic  energy  of  the 
smallest  grain  at  the  tapered  end  are  developed  for  various  chain  lengths  and  tapering  factors.  These 
details  are  then  compared  to  congruent  diagrams  obtained  via  extensive  dynamical  simulations. 
Our  figures  indicate  that  the  ratios  of  the  kinetic  energies  of  the  smallest  to  largest  grains  possess  a 
gaussian  dependence  on  tapering  and  an  exponential  decay  when  the  number  of  grains  increases.  The 
conclusions  are  independent  of  system  size,  thus  being  applicable  to  tapered  alignments  of  micron¬ 
sized  spheres  as  well  as  those  that  are  macroscopic  and  more  easily  realizable  in  the  laboratory.  The 
results  demonstrate  the  capability  of  these  chains  to  thermalize  propagating  impulses  and  thereby 
act  as  potential  shock  absorbing  devices.  While  inertial  mismatches  in  these  granular  chains  lead  to 
remarkable  energy  absorption,  short  chains  are  found  to  be  limited  in  that  regard.  A  second  granular 
system  is  therefore  proposed  and  investigated  which  greatly  improves  performance  for  any  size  chain. 
These  new  systems  feature  surprisingly  complicated  dynamics  and  are  inadequately  represented  by 
a  hard-sphere  approximation.  Additionally,  such  systems  have  shock  absorption  capacities  that  vary 
as  a  function  of  position  along  the  chain,  enabling  customized  shock  absorbers. 

Additional  studies  investigate  energy  partitioning  and  fluctuations  are  investigated.  Approximate 
power  laws  are  developed  which  fit  the  decay  of  average  fluctuations  as  the  size  of  the  system 
increases.  Advanced  simulations  of  tapered  chains  utilizing  the  modern  hyclrocode,  ALEGRA  is 
introduced.  These  simulations  incorporate  elastic-plastic  equation  of  state  and  behavior  allowing 


us  to  probe  very  large  loading  of  tapered  chains.  This  leads  into  the  discussion  of  our  continuing 
work  beyond  this  dissertation  including  the  design  of  a  shock  absorbing  panel.  Historical  context  is 
provided  which  has  lead  researchers  to  begin  looking  seriously  at  these  alluring  properties  of  granular 
or  discretized  systems. 


x 


List  of  Figures 


1.1  Hertz,  a:5/2,  Hookean-like,  x2,  and  hard  core  potentials  where  n  =  5, 10 .  11 

1.2  Phase  space  diagrams  for  equation  (1.8) .  15 

1.3  Single  particle  system  between  fixed  compressible  walls  under  harmonic  (n  =  2.0)  and 

anharmonic  (n  =  2.5)  potentials  (a)  normalized  kinetic  energy,  (b)  phase  space,  (c) 
force  versus  displacement .  17 

1.4  Half-periods  plotted  as  a  function  of  v-i  for  various  values  of  n.  Each  superimposed 
fit  has  the  form  given  by  T/2  =  A{n)v~x^  (dashed).  For  each  n  the  values  of  A,x 

are  plotted  in  Figure  1.5 .  18 

1.5  Coefficients  x,A  from  Figure  1.4  plotted  vs.  n  and  fitted  with  appropriate  functions 

in  support  of  generalized  expression  for  T/2 .  18 

1.6  Velocity  and  position  phase  diagrams  for  the  two  particle  system.  Dynamics  are 

shown  for  the  first  120  qis  with  corresponding  points  of  interest  identified.  Note  that 
maximum  velocities  do  not  necessarily  correspond  to  crossing  the  t  =  0  positions.  .  .  20 

1.7  Cartoon  depicting  the  evolving  two  particle  system.  Refer  to  figure  1.6  for  locations 

in  phase  space .  21 

1.8  Overlap,  <5,  as  a  function  of  time  for  the  particle-particle  interface  in  the  N  =  2 

system.  S  is  represented  as  a  percentage  of  the  two  spheres  combined  radii  since  it  is 
normalized  by  10  mm  (2 Rt) .  22 

2.1  Simple  tapered  chain:  N  =  10,  qs  =  8%,  L  =  70.7  mm,  ri  =  5  mm.  The  rightmost 

particle  is  grain  i,  its  nearest  neighbor  is,  i  +  1,  etc .  24 

2.2  Normalized  Ek(N,  q,  El)  parameter  space  for  the  STC  hard-sphere  approximation. 

N  varies  from  3  to  20  and  q  from  0%  to  10% .  27 

2.3  Normalized  Ek{N ,  q)  parameter  space  for  the  STC  hard-sphere  approximation  where 

initial  velocity  is  supplied  to  the  smaller  end  of  the  chain .  28 

2.4  KE  as  a  function  of  time  for  the  smallest  particle  in  a  chain  with  N  =  15,  u>  =  0.5. 

Each  plot  is  evaluated  for  a  different  tapering  with  initial  kinetic  energy  as  0.0838  J.  29 

2.5  KE  as  a  function  of  time  for  the  smallest  particle  in  a  chain  with  N  =  20,  ui  =  0.5. 

Each  plot  is  evaluated  for  a  different  tapering  with  initial  kinetic  energy  as  0.0838  J.  30 


2.6  Absolute  positions,  velocities,  and  kinetic  energies  for  particles  11-15.  Particle  15  is 

the  last  grain  and  is  in  contact  with  the  boundary.  The  subplots  (a-c)  represent  the 
monodisperse  chain  where  each  grain  has  r  =  5  mm,  while  (d-f)  corresponds  to  a 
chain  with  a  tapering  of  10%  so  that  particles  11-15  have  radii  1.74,  1.57,  1.41,  1.27, 
and  1.14  mm,  respectively.  These  subplots  illustrate  when  the  grains  in  question  first 
receive  the  incident  impulse.  Note  the  earlier  time  of  arrival  for  q  =  0.1  since  velocities 
are  higher  (up  to  a  factor  of  3  in  the  plots).  For  panels  (d-f)  the  dynamics  have  been 
extended  beyond  initial  incidence  and  reflection .  31 

2.7  Each  of  the  plots  depict  force  as  a  function  of  time  for  the  tail  and  head  particles  in  a 
chain  with  N  =  15,  u>  =  0.5.  The  value  of  the  tapering  increases  for  successive  plots 
and  is  denoted  by  the  label  q.  Negative  forces  imply  acceleration  towards  the  smaller 


end  of  the  chain .  32 

2.8  Numerical  solution  of  Ek{N,  q)  parameter  space  for  constant  u> .  36 


xi 


2.9  Numerical  solution  of  F(N,q)  parameter  space  for  constant  u> .  37 

2.10  Difference  plot  of  the  lossless  KE  parameter  spaces  for  the  hard-sphere  approximation 

(Figure  2.2a)  and  numerical  solution  (Figure  2.8a).  Note  that  the  azimuthal  view  has 
been  rotated  180  degrees  for  clarity .  37 

2.11  Energy  partitioning  for  STCs  where  N  =  {3, 8, 14, 20}  represents  the  number  of 

spheres  and  qs  =  {0.0,0.05,0.1}  is  the  tapering.  Noisy  plots  resembling  panel  l  in¬ 
dicate  shock  absorbing  systems  while  those  similar  to  panel  j  transmits  an  impulse 
as  a  solitary  wave  —  essentially  without  loss.  Panels  a-c  represent  efficient  energy 
conversion  systems .  38 

2.12  Normalized  KE  for  N  =  3,  qs  =  0.05.  Ei  represents  the  head  of  the  chain  and  input; 

E2  is  the  central  bead;  and  E$  is  the  tail.  Note  the  periodicity  of  the  system  resulting 

in  a  recurrence  time  of  just  under  80  /is .  39 

2.13  Instantaneous  kinetic  energy  per  grain  for  N  =  15,  t  =  52/us  and  selected  tapering.  .  40 

2.14  Normalized  kinetic  energy  landscape.  Vertical  axes  represents  the  particle  id  while 

the  horizontal  axis  is  time  out  to  1000  /is .  49 

2.15  Velocity  distribution  with  superimposed  Gaussian  fit  for  each  particle  in  a  STC  with 

N  =  4,  q  =  0.  Particle  1  represents  the  end  of  the  chain  and  particle  4  represents  the 
head  where  the  initial  impulse  was  applied .  50 

2.16  Velocity  distribution  with  superimposed  Gaussian  fit  for  each  particle  in  a  STC  with 
N  =  20,  q  =  0.  Particle  1  represents  the  end  of  the  chain  and  particle  20  represents 

the  head  where  the  initial  impulse  was  applied .  51 

2.17  Average  kinetic  energy  of  the  system  as  a  function  of  n  and  N .  52 

2.18  Average  kinetic  energy  of  the  system  determined  from  simulation  (circles)  and  theo¬ 
retical  expectation  by  the  Virial  Theorem  (dashed) .  52 

2.19  Partition  of  energy  by  grain  number  for  several  STCs  where  N  =  {3,8,14,20}  and 
q  =  {0, 0.05, 0.1}  at  t  =  100 /is.  The  red  line  indicates  the  theoretical  expectation  that 

<  KEn  >=  "/("+2) .  53 

2.20  Partition  of  energy  by  grain  number  for  several  STCs  where  N  =  {3,  8, 14,  20}  and 
q  =  {0, 0.05, 0.1}  at  t  =  500 /is.  The  red  line  indicates  the  theoretical  expectation  that 

<  KEn  >=  "/("+2) .  53 

2.21  Partition  of  energy  by  grain  number  for  several  STCs  where  N  =  {3,  8, 14,  20}  and 

q  =  {0,0.05,0.1}  at  t  =  1000 /is.  The  red  line  indicates  the  theoretical  expectation 
that  <  KEn  >=  "/("+2) .  54 

2.22  Fluctuations  as  a  function  of  time  and  n  for  N  =  100,  q  =  0 .  55 

2.23  (a)  Variation  of  <  F  >  with  N  for  various  n  and  fitting  function,  (b)  Variation 
(scatter)  of  k  with  n.  (c)  Variation  of  A  with  n  and  corresponding  fitting  function. 

(d)  Variation  of  B  with  n  and  corresponding  fitting  function .  56 

2.24  <  F  >  as  a  function  of  n  for  various  N  determined  from  simulation  (circles).  Derived 

fitting  function  <  F  >=<  F(n,  N )  >  is  superimposed  (dashed)  in  each  panel,  n  =  2 
data  is  excluded  as  it  doesn’t  fit  the  function .  57 

2.25  \vrnax\  as  a  function  of  time  for  N  =  16, 18,  20, 50 .  58 

2.26  The  Simulated  and  Modeled  STC .  59 

3.1  Several  example  DTCs  created  by  varying  /  and  qd  for  constant  N  =  13 .  61 

3.2  Normalized  kinetic  energy  surfaces,  Ek  =  Kout / Kin,  for  the  decorated  chain  under 
the  hard  sphere  approximation  as  functions  of  the  number  of  spheres,  N,  fractional 

size  of  interstitial  sphere,  /,  and  tapering,  qd .  65 

3.3  Numerically-produced  normalized  kinetic  energy  surfaces,  Ek  =  Kout / Kin,  for  the 

decorated  chain  as  functions  of  the  number  of  spheres,  N,  fractional  size  of  interstitial 
sphere,  /,  and  tapering,  qd .  Several  sample  chains  are  identified  in  panels  (d-i).  ...  67 

3.4  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  52/rs 

where  qd  =  {0,  0.05, 0.1}  and  /  =  {1.0, 0.7, 0.3} .  68 

4.1  Computational  domain  and  problem  setup  for  the  ALEGRA  simulation .  70 

xii 


4.2  Evolution  of  \crxx\.  Color  scale  is  logarithmic .  71 

4.3  Basic  3D  simulation  of  a  TC.  Note  that  the  bulk  behavior  of  the  spheres  is  consistent 

with  expectation  even  though  dynamics  take  place  on  a  elemental  or  nodal  basis.  A 
small  amount  of  plastic  flow  is  visible .  72 

4.4  Velocity  of  the  head  particle  as  a  function  of  time  for  both  EOM  and  ALEGRA 

simulations  where  N  =  5,  q  =  0.07  and  restitutive  losses  are  neglected .  73 

4.5  Time  elapsed  sequence  where  the  initial  velocity  of  the  left-most  grain  has  been  in¬ 
creased  to  500  m/s  (red).  The  remaining  spheres  are  initially  stationary  (blue).  ...  75 


5.1  Schematic  of  what  a  tapered  chain  armor  panel  might  look  like .  78 

D.l  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  10/us 

where  qd  =  {0,0.05,0.1}  and  /  =  {1.0, 0.7, 0.3} .  102 

D.2  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  20/us 

where  qd  =  {0,0.05,0.1}  and  /  =  {1.0, 0.7, 0.3} .  103 

D.3  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  30/xs 

where  qd  =  {0,  0.05, 0.1}  and  /  =  {1.0, 0.7, 0.3} .  104 

D.4  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  40/us 

where  qd  =  {0,  0.05, 0.1}  and  /  =  {1.0, 0.7, 0.3} .  105 

D.5  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  50/us 

where  qd  =  {0,0.05,0.1}  and  /  =  {1.0, 0.7, 0.3} .  106 

D.6  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  60/us 

where  qd  =  {0,0.05,0.1}  and  /  =  {1.0, 0.7, 0.3} .  107 

D.7  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  70/us 

where  qd  =  {0,0.05,0.1}  and  /  =  {1.0, 0.7, 0.3} .  108 

D.8  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  80/us 

where  qd  =  {0, 0.05, 0.1}  and  /  =  {1.0, 0.7, 0.3} .  109 

D.9  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  90/us 

where  qd  =  {0,  0.05, 0.1}  and  /  =  {1.0, 0.7, 0.3} .  110 

D.10  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  100/rs 

where  qd  =  {0, 0.05, 0.1}  and  /  =  {1.0, 0.7, 0.3} .  Ill 


xiii 


Chapter  1 

Introduction 


Most  of  everyday  life  is  nonlinear,  and  the  principle  of  superposition  fails  spectacularly.  If  you 
listen  to  your  two  favorite  songs  at  the  same  time,  you  wont  get  double  the  pleasure!  -  Steven  H. 
Strogatz 108 


1.1  Statement  of  the  Problem 

Most  everyone  is  familiar  with  or  has  seen  a  Newton’s  cradle.  This  popular  physics  apparatus  and 
desk  mainstay  marvelously  demonstrates  the  conservation  laws  of  energy  and  momentum.  A  further 
attribute  of  the  system  pointed  out  by  Hermann,  et.  al. 42-44  is  the  requirement  of  dispersionless 
propagation  when  there  are  more  than  two  spheres.  Thus,  the  spheres  must  be  identical.  What 
most  people  may  not  be  aware  of  however,  is  that  by  tapering  the  size  of  the  spheres  the  system 
becomes  a  shock  absorber. 

This  has  broad  implications  since  shock  mitigation  is  one  area  that  will  always  be  receptive  to 
improvements  in  the  the  state  of  the  art.  It  encapsulates  several  important  applications  of  military, 
commercial,  and  industrial  interest,  such  as  blast-proofing,  vibrational  or  sound  suppression,  and 
noise  filtering.  Traditional  methods  of  dealing  with  undesirable  transients,  such  as  that  from  ballistic 
shock,  include  metal  foams  and  honeycombs33,34.  When  honeycombs  are  extruded,  one  obtains  a 
linear  celluar  alloy41,  which  has  demonstrated  improved  energy  absorption  capabilities32.  Another 
approach  to  dismiss  transients  is  through  the  use  of  functionally  graded  materials  (FGM) 14  where 
one  can  introduce  impedance  mismatches  gradually  or  discontinuously. 

The  general  problem  under  investigation  therefore  is  understanding  the  energy  mitigation  process 
behind  one-dimensional  alignments  of  tapered,  metal  spheres  in  contact.  In  particular,  this  has  led 
to  the  investigation  of  two  families  of  systems  which  are  broadly  referred  to  as  tapered  chains  (TC). 
These  alignments  barely  touch  in  their  initial  configuration  and  are  not  under  any  precompression. 
In  addition,  they  are  considered  to  be  maintained  between  two  fixed  —  but  compressible  —  walls. 
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When  an  impulse  is  applied  to  one  end  of  the  system,  it  propagates  via  interparticle  contacts  and  can 
be  dissipated  or  maintained  depending  on  the  number  of  spheres,  the  ratio  of  adjacent  particle  sizes 
and  energy  losses.  Energy  coming  out  of  the  system  (or  anywhere  therein)  can  be  easily  measured 
in  computer  simulations  thus  allowing  one  to  quantify  the  energy  absorption  capability  as  a  function 
of  the  system  parameters.  In  this  work,  components  of  the  system  may  be  referred  to  as  grains, 
particles,  beads,  or  spheres. 

It  will  be  seen  that  such  seemingly  simple  systems  exhibit  very  interesting  nonlinear  dynamics. 
Indeed,  the  potential  energy  is  strongly  nonlinear  and  the  equations  of  motion  (EOM)  have  no 
known  analytic  solution.  Even  if  there  were  a  solution  in  closed  form,  it’s  not  clear  that  one  could 
divine  any  physical  meaning  or  intuition  from  it  due  to  its  presumed  complexity.  As  a  result,  and 
as  is  typical  for  nonlinear  problems,  one  resorts  to  numerical  methods  to  solve  for  the  particles’ 
position,  velocity,  and  acceleration  as  a  function  of  time.  At  this  point,  it  is  instructive  to  provide  a 
history  of  the  developments  which  have  led  us  to  this  study — an  amalgam  of  granular  media,  contact 
mechanics,  and  one-dimensional,  discrete  systems. 

1.2  Background 

The  propagation  of  mechanical  impulses  along  a  chain  of  spheres  has  become  an  increasingly  strong 
research  area  during  the  past  decade.  Tapered  chain  systems  are  customized  ID  constructions 
of,  fundamentally,  granular  media  where  the  contact  nature  is  quite  intriguing  and  critical  to  the 
dynamics.  Granular  media  is  usually  considered  to  be  2D  and  3D  entities.  A  brief  review  of  some 
observed  phenomenology  (section  1.2.1)  is  useful  since  similar  behavior  may  be  seen  and  prove  useful 
in  the  analysis  of  ID  systems.  Section  1.2.2  provides  examples  of  the  computational  and  theoretical 
studies  of  granular  chains,  while  section  1.2.3  covers  the  experimental  endeavors. 

1.2.1  Granular  Media 

Discrete  or  granular  media29, 48,49,51,64,91  consists  of  particles  that  can  range  in  size  from  microm¬ 
eters  to  meters  and  number  in  quantity  from  several  to  the  uncountable.  The  most  well-known 
constituents  of  this  group  are  sands  and  powders  and  are  utilized  across  many  disciplines.  This 
media  is  particularly  intriguing  in  that  among  other  things  it  has  properties  of  liquids  and  solids 
and  calls  have  been  made  to  add  it  as  a  fundamental  state  of  matter  29 .  It  exhibits  alluring  behavior 
is  seen  in  the  phenomena  of  force  networks,  avalanches,  and  jamming.  Mathematically,  granular 
media  can  be  shown  to  be  similar  to  nonlinear  rods120. 
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Granular  media  (in  any  number  of  dimensions)  has  been  intensely  investigated  for  understanding 
wave  propagation1.  Much  of  this  has  been  in  regard  to  understanding  sound  propagation.  It  has 
become  such  a  curiosity  since,  in  2D  or  3D,  not  all  grains  are  under  contact  so  various  contact  or 
force  networks  are  formed.  That  is,  only  particles  in  contact  admit  mechanical  energy  propagation 
and  the  identified  network  visually  appears  dendritic.  However,  when  one  slightly  shakes  the  gran¬ 
ular  bed,  the  network  is  completely  changed.  This  is  important  for  improving  imaging  of  buried 
objects,  such  as  munitions,  in  various  soils.  Sinkovits  and  Sen101,102  presented  simulations  on  the 
vertical  propagation  of  weak  and  strong  impulses  in  deep  gravitationally  compacted  columns.  In 
such  asymmetrically  loaded  columns,  the  sound  velocity,  c,  increases  with  depth  as  c  oc  z1/6 ,  where 
z  is  the  depth.  Sen101  extended  the  study  to  include  voids  and  mass  impurities.  The  inclusion  of 
the  latter  caused  backscatter  at  the  defects  even  though  the  chains  were  monodisperse  (identically- 
sized).  This  then  prompted  further  study  by  Sen,  et.  al.9S  to  report  it  as  a  possible  mechanism  for 
the  detection  of  buried  impurities,  such  as  mines,  bones,  lost  treasure,  etc.  The  specific  nature  of 
the  propagation  59  and  backscattering60  was  investigated  more  thoroughly  by  Manciu  who  correlated 
the  results  into  a  dissertation58  on  nonlinear  acoustics  in  granular  media. 

It  is  intriguing  and  suggestive  that,  conceptually,  granular  media  can  be  envisioned  as  the  inverse 
of  a  porous  material.  In  a  porous  material,  there  are  gaseous  (air)  voids  within  a  solid  matrix.  For 
dry  granular  media,  the  reverse  is  true:  solid  (grains)  voids  exist  within  a  gaseous  matrix.  This  is 
an  interesting  comparison  to  draw  because  the  considerable  amount  of  work  required  to  collapse  the 
voids  in  a  porous  material  makes  it  a  favorable  technology  against  ballistic  shock32. 

1.2.2  Granular  Chains:  Computational  and  Theoretical  studies 

Since  the  Hertz  potential  is  strongly  nonlinear,  a  review  of  nonlinear  methods  and  waves  was  in 
order.  It  was  believed  that  such  resources  would  prove  valuable  in  the  analysis  and  approxi¬ 
mation  of  TCs.  While  there  are  an  inordinate  number  of  books  available  on  nonlinear  dynam¬ 
ics  39’40,67’68’73’74,105>106i108,ii4^  they  generally  focus  on  perturbation-type  solutions  and  analysis  of 
the  resulting  nonlinear  and  chaotic  waves9, 70’106,112,118. 

In  most  cases,  undamped,  damped  and  forced  oscillators,  such  as  the  Van  der  Pol  and  Duffing 
oscillators,  are  examined  in  detail  along  with  the  their  behaviors  in  phase  space.  The  behavior  of 
nonlinear  waves  is  introduced  through  the  simplest  forms  of  dispersive  and  diffusive  equations9,73: 
the  Korteweg  deVries  (KdV)  equation52,122  and  Burgers  equation15,  respectively.  The  KdV  equation 
for  example  admits  solitary  wave  solutions. 
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A  shortcoming  of  these  works  is  that  none  of  them  address  fractional  power  exponents.  Mickens67 
and  Gottlieb38  have  the  most  relevant  analytic  work  on  oscillators  with  fractional  power  nonlinear¬ 
ities.  Their  approach  uses  the  so-called  Harmonic  Balance  technique66  for  periodic  systems  which 
matches  the  coefficients  in  a  truncated  Fourier  series  to  determine  the  frequency  of  oscillation.  A 
benefit  of  the  method  is  that  it  does  not  require  a  perturbation  parameter  so  it  is  applicable  to 
strongly  nonlinear  systems.  Harmonic  Balance  however  seems  to  be  limited  to  specific  problems. 

Solitary  waves  (SW)  are  intriguing  because  they  propagate  over  very  large  distances  with  negligi¬ 
ble  dissipation86.  The  first  known  observation  was  reported  in  1844  by  Russell92  who  noticed  water 
waves  in  a  shallow  canal  propagate  without  loss  over  several  miles.  Ultimately  they  were  broken 
apart  by  branching  of  the  channel.  An  exhaustive  historical  account  of  the  theoretical  development 
of  SWs  is  expounded  by  Sander93.  It  should  be  noted  that  when  SWs  collide  and  emerge  unchanged 
they  are  referred  to  as  solitons  because  of  the  similarity  to  particles83. 

SWs  were  of  particular  interest  in  past  studies  since,  for  monodisperse  chains,  they  are  the 
mode  of  energy  transport77,95.  Sen  and  Pfannes  demonstrated  that  it  takes  about  15  grains  for 
the  SW  to  form  in  such  a  chain.  Considerable  effort  has  been  spent  investigating  the  phenomena 
numerically  and  experimentally.  Nesterenko  first  arrived  at  an  analytic  solution  of  a  precompressed, 
monodispersed  system  using  anharmonic  and  long  wavelength  approximations  of  the  EOM75.  The 
equations  were  then  reduced  to  a  form  similar  to  the  KdV  equation52  whereby  SWs35  were  observed. 

Lattice  dynamics  are  dependent  on  the  type  of  interaction  potential  chosen.  The  ID  (atomic) 
lattice  and  coupled  oscillators  immediately  stand  out  where  Morse7,  exponential114,  and  other 
nonlinear  (pertubative)  potentials  have  been  used.  A  common  observable  in  most  cases  is  the 
existence  of  solitary  waves.  In  point,  Toda114  obtained  exact  (soliton)  solutions  in  a  one-dimensional 
lattice  with  exponential  interactions  in  closed  form. 

Additionally,  one  of  the  more  curious  behaviors  in  granular  chains  is  that  of  inelastic  collapse63. 
This  occurs  when  the  separation  between  particles  drops  to  zero  due  to  a  competition  between  the 
number  of  particles  in  the  system  and  energy  losses.  Since  the  velocity  of  the  particles  involved  in 
the  collapse  decay,  this  would  be  an  ideal  phenomenon  to  exploit  for  shock  mitigation  purposes.  We 
have  not  yet  investigated  this. 

In  the  following  subsections,  specific  developments  of  FPU  and  mixed  harmonic-anharmonic 
potentials  and  the  overlap  potential  are  discussed.  This  follows  with  some  review  of  work  focused 
on  equipartition  and  equilibrium  in  such  systems. 
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FPU  and  mixed  harmonic-anharmonic  potentials 


Fermi,  Pasta,  and  Ulam30  (FPU)  are  credited  with  the  first  computational  study  on  energy  sharing 
between  modes  where  quadratic,  cubic,  and  broken  linear  potentials  were  used.  As  a  side  note,  the 
FPU  ID  lattice  has  become  a  testbed  for  looking  at  various  physics  such  as  impulse  propagation, 
defining  temperature,  and  energy  sharing  among  modes.  An  enormous  body  of  literature  is  available 
on  the  subject.  They  noticed  that  upon  the  inclusion  of  such  nonlinearities,  energy  that  was  initially 
fed  to  the  lower  mode  became  mixed  among  several  lower  order  modes,  but  then  returned  to  the 
lowest  mode.  Such  recurrence  can  indicate  hidden  periodicities.  These  results  are  in  contrast  to  a 
purely  linear  problem.  For  example,  if  one  were  to  excite  the  lowest  mode  only  in  a  general  solution 
to  the  (linear)  wave  equation,  the  energy  would  stay  in  the  fundamental  mode  for  all  time78. 

The  approach  to  equilibrium  in  discrete  systems  has  been  a  topic  of  considerable  interest  for  quite 
some  time  and  the  open  literature  is  replete  with  examples.  In  general  one  wishes  to  understand 
whether  the  system  reaches  equilibrium  by  how  energy  is  shared  among  the  sites  in  a  ID  lattice 
and  how  long  it  takes  to  do  so.  Typically  the  systems  studied  have  a  Hamiltonian  of  the  form, 
H  =  Hq  +  eHi  where  Hq  represents  the  linear  portion,  Hi  represents  nonlinear  (sin(;r),  x 3,  x4, 
etc.)  terms  and  e  is  a  tuning  parameter.  In  this  sense,  many  of  the  perturbation  techniques  in 
previous  references  may  be  applied.  An  additional  benefit  is  that  one  can  monitor  how  energy  is 
shared  between  normal  modes  from  the  linear  system  and  the  effects  of  adding  a  small  nonlinear 
perturbation. 

Fermi  addressed  the  question  of  energy  equipartition 30  and  Toda114  points  out  lucidly, 

Fermi  did  some  work  on  the  ergodic  problem  when  he  was  young,  and  when  electronic 
computers  were  developed  he  came  back  to  this  as  one  of  the  problems  computers  might 
solve.  He  thought  that  if  one  added  a  nonlinear  term  to  the  force  between  particles  in 
a  one-dimensional  lattice,  energy  would  flow  from  mode  to  mode  eventually  leading  the 
system  to  a  statistical  equilibrium  state  where  the  energy  is  shared  equally  among  the 
modes. 

This  investigation  sparked  a  broad  literature  base24,85  and  a  description  of  some  of  this  work 
follows.  Ford  and  Waters  have  addressed  the  issue  of  equipartition  of  FPU  systems  in  several 
exhaustive  papers.  Ford31  used  a  perturbation  technique  to  show  that  the  chain  originally  used 
by  FPU  didn’t  equipartition  energy  because  of  an  inappropriate  choice  of  frequencies.  He  argues 
that  for  weakly  coupled  oscillator  system  —  with  linear  and  nonlinear  terms  —  one  should  expect 
(internal)  resonance  and  energy  sharing  when  the  sum  of  harmonic  frequencies  is  approximately 
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zero. 


Tobolsky  et.  al. 113  observed  equipartition  for  several  lattice  models.  They  commented  however 
that  the  probability  of  a  given  deviation  of  kinetic  and  potential  energies  from  one-half  approaches 
a  limit  that  decreases  with  increasing  N.  For  a  deviation,  x,  with  probability  of  1%  being  exceeded, 
I  measure  this  decay  (based  on  their  data)  as  x  =  3/(4 VN),  where  N  is  the  number  of  particles. 

Boccheri  et.  al. 10  numerically  investigate  a  one  dimensional  chain  of  particles  with  nearest  neigh¬ 
bor  interactions  using  a  Lennard-Jones  potential.  They  conclude  that  when  the  particle  vibration 
energy  exceeds  a  few  percent  of  the  depth  of  the  potential  well,  equipartition  of  energy  among  the 
normal  modes  is  satisfied  in  the  time  average.  This  is  in  contrast  to  the  Toda  lattice114  where 
exponential  interactions  were  not  compatible  with  equipartition  of  energy  in  the  time  average. 

Rosas  and  Lindenberg89  studied  pulse  propagation  in  FPU  lattices  (where  only  quartic  nonlin¬ 
earities  were  considered) .  They  varied  the  input  velocity  governing  whether  harmonic  or  anharmonic 
terms  dominate  (since  the  velocity  of  nonlinear  waves  is  amplitude  dependent).  They  concluded  that 
the  pulse  width  is  not  an  appropriate  measure  of  the  way  a  pulse  spreads  in  a  purely  anharmonic 
lattice  but  rather  it  measures  the  span  over  which  a  series  of  decreasing  velocity  pulses  exist.  This 
was  an  extension  of  an  earlier,  more  comprehensive,  study  by  Sarmiento  et.  al . 94.  They  measured 
pulse  propagation  in  isolated  FPU  chains  and  also  those  coupled  to  heat  baths  of  zero  and  finite 
temperatures  as  well  as  2d  isolated  arrays. 

Finally,  one  of  the  more  recent  discoveries  has  been  the  observation  of  “breathers”  -  -periodic  and 
localized  nonlinear  lattice  excitations81.  These  spontaneous  and  long-lived  localizations  of  energy 
result  from  competition  between  nonlinearity  and  space  discreteness81.  Another  typical  requirement 
is  for  the  array  to  be  cooled  at  the  boundaries.  Essentially,  all  energy  dissipation  takes  place  at  the 
end  points  and  is  equivalent  to  the  physical  situation  of  a  surface  cooling  much  faster  than  the  bulk 
of  the  medium.  Breathers  emerge  from  such  configurations  where  the  lattice  is  initially  thermalized 
and  is  an  evolved  form  of  self-excited  oscillations — another  nonlinear  behavior.  In  this  report, 
such  behavior  has  not  been  observed  because  free-encl  BCs  nor  boundary-only  cooling  have  been 
considered.  Further  conclusions  of  Piazza,  et.  al.  were  that  breather  mobility  is  impacted  if  space 
discreteness  (i.e.  weak  nearest  neighbor  interaction  or  restoration  potential)  is  too  large. 

Overlap  potential 

The  study  of  granular  columns  is  dominated  by  the  Hertz  potential  —  see  equation  (1.3).  This  is 
a  special  case  of  the  more  generalized  overlap  potential,  6n.  The  overlap  potential  is  quite  distinct 
from  FPU  and  other  lattice  potentials.  First,  they  lack  a  harmonic  term  —  that  is,  they  are  strongly 
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nonlinear.  Thus,  nonlinear  perturbation  methods  fail  when  looking  for  a  solution.  Second,  and  most 
important,  is  that  they  lack  a  restoring  term.  So  even  if  the  overlap  potential  is  considered  with  a 
“harmonic”  exponent,  <5  ,  the  chain  is  very  dissimilar  from  a  simple  harmonic  system. 

Nesterenko  is  considered  to  be  the  pioneer  of  theoretical,  computational,  and  experimental  studies 
of  SW  in  Hertzian  chains.  His  work  is  found  compiled  in  chapter  one  of  his  book 76  which  includes 
translations  of  his  Russian  publications.  He  was  the  first  to  observe  solitons  in  Hertzian  chains  '5 
and  also  dubbed  such  granular  chains  a  “sonic  vacuum,”  since  particles  can  lose  contact  when 
precompression  in  the  chain  is  not  considered. 

The  phenomenology  of  SWs  in  granular  chains  was  addressed  specifically  in  several  computa¬ 
tional  papers  by  Sen  and  Manciu 61,95 .  For  long,  uncompressed  chains — where  boundary  effects 
are  ignored — they  observed  that  SWs  are  the  mode  of  energy  transport  when  n  >  2.  They  used 
a  combination  of  computational  and  analytic  arguments  to  guess  a  form  of  the  SW  as,  4>n(z)  = 
—  {A/ 2)  tanh[/„(z)/2],  where  the  fn(z )  represents  a  series  expansion  in  2  and  requires  calculation  of 
the  coefficients  for  each  term  in  the  series.  Their  solution  agrees  quite  well  with  simulations  when 
only  the  first  three  coefficients  are  solved.  Thus  they  were  able  to  obtain  displacement,  velocity,  and 
acceleration  functions  for  the  SW.  This  approach  differs  from  that  of  Nesterenko’s  long  wave,  or  con¬ 
tinuum,  approximation  in  that  zero  precompression  is  considered  and  the  EOM  are  not  linearized. 
Their  subsequent  publication  investigates  how  propagation  of  the  SW  is  affected  by  dissipation  and 
disorder  for  arbitrary  power-law  repulsive  potentials,  6n.  They  conclude  that  even  in  the  presence  of 
disorder,  a  SW  propagates — maintaining  its  width  while  the  amplitude  decreases  exponentially  with 
distance.  The  attenuation  is  common  whether  energy  loss  is  caused  by  restitution  or  velocity-based 
friction,  or  from  multiple  backscattering  events  due  to  randomly-sized  particles  in  the  chain.  The 
attenuation  commonality  thus  enables  one  to  choose  an  energy  loss  mechanism  that  is  convenient. 
An  additional  note  of  interest  is  that  as  the  SW  propagates  in  the  disordered  medium,  part  of  the 
energy  remains  behind  the  leading  edge  in  the  form  of  noise  which  doesn’t  interfere  with  the  SW. 
Since  wave  velocity  is  a  function  of  amplitude  and  noise  has  a  much  smaller  amplitude,  the  SW 
outpaces  noise. 

In  a  subsequent  study,  Sen  et.  al. 97  inquired  as  to  whether  it  was  possible  to  convert  an  impulse 
into  thermal  energy.  Based  on  supporting  computational  studies  they  concluded  that  when  particle 
radii  taper  to  smaller  sizes,  the  traveling  SW  must  get  “squeezed”  into  a  smaller  size.  The  SW 
consequently  loses  its  reflective  symmetry  and  is  destroyed.  More  specifically,  the  SW  is  rapidly 
attenuated  and  undergoes  dispersion  due  to  progressively  smaller  and  faster  particle  masses.  That 
the  energy  transport  mechanism  was  found  to  be  disabled  by  tapered  chains  led  to  many  exciting 


7 


questions  and  became  the  foundation  for  this  dissertation.  Curiously,  Poschel  and  Brilliantov 82 
found  that  by  using  a  constant  coefficient  of  restitution,  one  obtains  optimum  energy  transmission 
if  the  mass  of  each  particle  is  prescribed  by  an  exponentially  decreasing  function. 

Wu121  uses  an  independent  collision  or  binary  collision  model  to  understand  wave  propagation 
in  uniform  (monodisperse)  and  tapered  chains.  In  such  a  system,  energy  is  always  confined  to  a 
single  particle.  This  model  is  analogous  to  hard-spheres  which  is  derived  in  sections  2.2  and  3.2. 

Rosas  and  coauthors  perform  a  series  of  relevant  studies  with  Hertzian  chains.  In  the  first  work87, 
the  authors  look  at  how  propagation  and  backscattering  between  two  granules  is  affected  by  the  type 
(kinetic  or  hydrostatic)  of  friction  and  its  magnitude.  Among  their  conclusions  is  that  friction  on 
the  second  granule  is  responsible  for  backscattering.  And  in  sum,  they  found  a  drastic  asymmetry 
between  backscattering  and  propagation  for  various  frictional  states.  This  work  is  followed  up88  by 
concluding  that  a  binary  collision  model  is  quantitatively  correct  for  (hard)  potentials  where  n  >  3. 
This  is  because  the  pulse  is  so  narrow  that  at  any  given  moment  the  energy  is  concentrated  in  just 
a  few  particles.  Indeed  they  reiterate  that  in  the  continuum  approximation  ( N  — >  oo),  the  pulse 
width  goes  as  a  =  \J (6 (n  —  2 ))2/n(n  —  1).  Subsequently,  the  authors  investigate  frictional  effects 
further  for  cylindrical  and  spherical  granules90.  An  interesting  result  for  the  cylindrical  case  is  that 
the  backscatter  velocity  for  finite  friction  is  greater  than  that  when  it  is  neglected.  They  also  find  an 
exponential  decay  to  the  velocity  of  backscatter  velocity  as  well  as  energy.  The  latter  is  in  agreement 
with  findings  in  this  report  (see  figure  2.26(a) — T(u>)  decays  exponentially) 

In  more  recent  work,  a  new  type  of  steady-state  behavior,  quasi-equilibrium99,100,  has  been 
observed.  The  findings  of  the  authors’  are  summarized  as  follows.  Without  a  linear,  harmonic  term 
in  the  potential,  grains  do  not  exhibit  simple  harmonic  motion  which  produces  sound  waves  — 
indeed  a  restoring  term  is  required.  This  is  the  “sonic  vacuum”  Nesterenko  referred  to.  As  such, 
information  or  energy  is  transmitted  by  groups  of  particles  rather  than  individual  grains.  When 
an  interaction  with  a  boundary  occurs,  energy  remains  nucleated  at  the  site,  some  of  which  goes 
into  recreating  an  attenuated  SW  while  the  rest  rattles  about,  creating  secondary  SW  (SSW) 57  of 
very  small  amplitude.  Since  SW  width  is  controlled  by  n,  the  system  transmits  energy  via  SW 
and  SSWs  such  that  (1)  large  and  small  SW  collisions  exchange  momenta  if  parallel114;  (2)  they 
undergo  breakdown  and  reconstruction  with  attenuation  by  antiparallel  collisions  or  interaction 
with  a  boundary59.  At  large  times,  it  was  found  that  the  system  tends  towards  a  Gaussian  velocity 
distribution.  Apparently  in  all  cases  studied,  the  process  of  reconstructing  a  SW  is  imperfect.  There 
is  attenuation  every  time  a  boundary  is  encountered,  the  remaining  energy  goes  into  making  SSW. 

Mohan  and  Sen  69  investigated  similar  questions  for  mass-spring  systems  with  quartic  interactions 


for  both  periodic  and  rigid,  perfectly  reflecting  boundary  conditions.  They  found  that  for  the  purely 
harmonic  mass-spring  chain,  energy  quickly  approaches  Eq/N  per  particle,  thus  equipartition  was 
achieved.  In  the  mixed  harmonic/anharmonic  case,  the  presence  of  harmonicity  tends  to  “wash¬ 
out”  nonlinear  effects.  Therefore  the  comparative  strengths  of  linear  and  nonlinear  terms  become 
an  important  factor  in  determining  the  relaxation  time. 


1.2.3  Granular  Chains:  Experimental  Work 

Some  of  the  earliest  experimental  work  focused  on  the  establishment  of  the  restitutive  coefficient 
for  various  materials  and  contact  time  of  impacts.  For  example,  Goldsmith36  illustrated  that  the 
duration  of  contact,  r,  for  two  identical  materials  is, 
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which  is  consistent  with  Hertz  theory. 

The  validity  of  the  Hertz  potential  has  been  discussed  in  several  papers21,75.  Bokor  and  Lev- 
enthall11  verified  Hertz  validity  outside  its  theoretical  regime — mainly  when  permanent  (inelastic) 
surface  deformations  occurred.  Experiments  have  also  been  performed  which  look  at  force  versus 
time  plots46  of  colliding  spheres  to  determine  the  power  of  the  restoring  term. 

Rossmanith  and  Shukla91  performed  the  first  study  on  photoelastic  investigations  of  granular 
media.  Their  work  spotlighted  dynamic  load  transfer  in  ID  vertical  and  increasingly  oblique  zig¬ 
zag  columns  of  the  photoelastic  material,  Homalite  100.  This  was  further  expanded  to  2D  beds  of 
roughly-surfaced  and  multi-sized  disks.  The  reported  isofringe  patterns  in  general  correspond  to 
higher  stress  accumulations.  In  addition,  their  report  represents  an  excellent  validation  exercise  for 
3D  hydrocode  simulations  (4). 

Much  experimental  work  has  been  carried  out  by  Nesterenko  and  collaborators.  One  interesting 
area  receiving  attention,  is  through  clever  adjustment  of  material  parameters.  In  this  case,  beads 
alternate  in  material  composition  rather  than  size.  This  has  been  investigated  experimentally  by 
Daraio23,  yet  I  am  unaware  of  numerical  studies. 

Coste  et.  al. 20  experimentally  validate  Nesterenko’s  work  where  both  zero  and  finite  precompres¬ 
sion  is  considered  through  force  measurements  on  a  sensor.  The  authors  stress  that  no  adjustable 
parameter  had  been  used  to  corroborate  their  results. 

Warr  and  Huntley11'  derive  from  experiment  an  exact  expression  for  the  rate  of  energy  input  into 
a  system  consisting  of  a  single  particle  vibrating  on  a  base  plate  in  one  dimension.  The  aggitation 
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coupled  with  the  rate  of  energy  dissipation  allowed  the  steady  state  of  the  system  to  be  determined 
energetically.  In  particular  they  found,  Ep  =  1.7mpV2 /(I  —  e),  where  V  is  the  base  plate  velocity, 
e  the  coefficient  of  restitution,  and  Ep  and  mp  are  the  particle  energy  and  mass,  respectively.  Ad¬ 
ditional  research56  investigates  the  phase  behavior — fluidization  and  condensation  (also  known  as 
clumping) — of  a  column  of  beads  under  similar  vibrating  plates. 

Nakagawa  et.  al.  '2  experimentally  confirm  impulse  dispersion  in  tapered  chains.  They  further 
provide  a  derived  measurement  for  the  coefficient  of  restitution  of  0.95  and  account  for  discrepancies 
in  the  computational  studies  by  arguing  that  some  of  the  energy  is  transferred  to  rotational  degrees 
of  freedom.  Experimental  work  has  also  focused  on  the  observation  of  SWs  and  their  behavior  at 
boundaries50,65’104.  With  relevance  to  chapter  3,  a  “decorated  tapered  chain”  prototype  has  very 
recently  been  constructed  for  empirically  testing  and  validating  its  shock  mitigation  capability2. 

1.3  Mathematical  Description  of  Problem 

We  define  TCs  (see  for  example  figure  2.1)  as  1-D  granular  arrays  of  elastic  spheres  that  touch  at  a 
single  point  in  their  initial  state  and  grow  to  a  disk  under  compression  in  the  plane  perpendicular 
to  the  figure.  The  chains  can  be  characterized  by  the  number  of  grains,  N,  the  successive  decrease 
in  size  of  the  grains  or  tapering,  q ,  and  restitutive  or  energy  losses,  u>. 

The  Hamiltonian  of  the  system  is  represented  as, 

#  =  ^Eto^+Ef(5m+ i)  (L2) 

i  i 

where  V(<5" i+1)  is  the  overlap  potential.  When  n  =  5/2,  V  is  referred  to  as  the  Hertz  potential  and 
is  described  in  section  1.3.1.  The  equations  of  motion  follow  in  section  1.3.2,  resititutive  losses  in 
section  1.3.3,  relevant  boundary  and  initial  conditions  in  section  1.3.4,  and  the  numeric  approach  to 
solving  the  N  equations  of  motion  in  section  1.3.5. 

1.3.1  The  Hertz  Potential 

The  contact  mechanics  between  adjacent  elastic  spheres  was  first  identified  by  Hertz45  circa  1882. 
Derivations  of  the  interaction  potential  were  later  performed  by  both  Landau53  and  Love55,  with 
a  simplified,  order  of  magnitude  approach  by  Leroy54.  The  results  were  that  for  spheres  under 
compression,  one  obtains  a  completely  repulsive  and  nonlinear  potential.  This  Hertz  potential  can 
be  written ,5,97  as, 
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Figure  1.1:  Hertz,  a;5/2,  Hookean-like,  x2 ,  and  hard  core  potentials  where  n  =  5, 10. 
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where 


(1.4) 


Here,  <5^+1  =  +  Ri+i  —  (^i+i  —  £j)  >  0,  represents  the  overlap  between  successive  grains  where 

Zj  is  their  position.  Additionally,  the  constant  a^i+i  has  been  defined  for  material  properties:  Ej,  the 
Young’s  modulus  and  <jj ,  the  Poisson  ratio;  and  radii,  R:r  Note  that  j  can  refer  to  either  particle 
i  or  i  +  1.  The  use  of  an  overlap  function  is  to  supplant  the  complicated  details  of  compression 
and  expansion.  Note,  that  the  nonlinearity  in  the  potential  is  completely  due  to  geometric  effects. 
Further,  the  theory  does  not  take  into  account  plasticity. 

For  our  particular  study,  all  materials  are  the  same  in  any  single  tapered  chain  so  that  equation 

(1.4)  reduces  to: 


(1.5) 


If  <  0  then  V  =  0  since  adjacent  grains  i  and  i  +  1  have  lost  contact.  It  may  be  noted 
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that  equation  (1.3)  describes  a  repulsive  potential  that  grows  faster  than  a  quadratic — or  Hookean- 
like — form  of  <52i+1.  The  phrase,  “Hookean-like”  is  used  to  stress  the  lack  of  a  restoring  term.  The 
Hertz  repulsion  is  hence  a  nonlinear  force.  More  specifically,  the  repulsion  is  softer  than  something 
harmonic  over  short  distances,  but  becomes  steeper  than  a  harmonic  form  with  increasing  compres¬ 
sion.  An  extreme  case  would  be  the  hard-sphere  potential  (n  — >  oo).  These  details  are  illustrated  in 
Figure  1.1  where  V{x)  is  plotted  for  various  values  of  n. 


1.3.2  Equations  of  Motion 

The  EOM  for  grain  rrii  at  position  Zj  is  constructed  from  equation  (1.3)  as 


"  _  5  r  .3/2  ^  x3/2  \ 

izi  —  2  ,i  p 


(1.6) 


where  the  dots  imply  differentiation  with  respect  to  time.  Recall  that  (pi+i  =  Ri  +  Ri+i  —  (2*+ 1  —  z-i) 
represents  the  overlap  of  successive  grains  where  Zj  is  the  position  of  a  grain.  Results  have  been 
obtained  for  a  large  selection  of  chains  consisting  entirely  of  TigALV  and  SiC  spheres.  Arbitrarily, 
we  have  chosen  to  use  TpARV  when  restitutive  losses  are  ignored  (ui  =  0),  and  SiC  otherwise.  The 
following  material  properties62  were  assumed  (where  D  is  defined  in  equation  (1.5)): 


Table  1.1:  Material  Properties 


Material 

p(mg/mm'3) 

D(mm2/kN) 

Occurence 

SiC 

3.2 

0.003266 

iO  0 

Ti6Al4V 

4.42 

0.01206 

u>  =  0 

Note  that  the  fundamental  units  of  length,  time,  mass,  and  force  in  our  simulations  are  the 
millimeter,  microsecond,  milligram,  and  kilonewton,  respectively.  This  is  to  minimize  potential 
rounding  errors  from  small  numbers. 

1.3.3  Restitutive  Losses 

Real  systems  have  various  modes  of  energy  dissipation — sliding,  rolling,  sound,  etc.  The  literature 
is  replete  with  methods  for  invoking  restitutive  losses5,6-19,  many  of  them  velocity  based109.  For 
simplicity,  we  have  utilized  the  method  of  Walton  and  Braun116,  where  the  coefficient  of  restitution, 
uj,  is  defined  as  Funioad/ Fioad  =  1  —  w.  Here  Funioad  represents  the  expansion  phase  of  the  contact 
event,  and  Fioad  is  the  compression  phase.  Values  of  u>  are  constant  (thus  hydrostatic)  throughout 
the  simulation  for  each  TC  and  are  chosen  as  0  <  u>  <  0.1.  Therefore  perfectly  elastic  collisions 
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correspond  to  u>  =  0.  The  majority  of  this  report  however  focuses  on  u>  =  0  since  it  establishes  an 
upper  limit  or  worst  case  scenario  in  terms  of  energy  mitigation. 

1.3.4  Boundary  and  Initial  Conditions 

In  our  model,  the  two  boundaries  consist  of  fixed,  compressible  walls.  This  is  equivalent  to  spheres 
of  infinite  radius.  As  a  result,  the  potential  is  adjusted  such  that  Rq,Rn+i  — >  oo.  At  an  interface 
with  the  boundary,  equation  (1.3)  becomes 

Vb  =  ^(R-z^,n)5/2, 

where  R  and  2  represent  a  particle  adjacent  to  the  boundary.  One  should  therefore  expect  behavior 
at  the  boundary  to  be  much  different  than  that  for  interior  particles. 

The  initial  conditions  have  been  based  on  a  delta  impulse  applied  to  a  bead  at  the  edge  (head) 
of  the  chain.  Specifically,  the  conditions  are 

Xi=N (t  =  0)  yf  0, 

=  0)  =  0. 

Unless  otherwise  stated,  the  Nth  particle  is  given  an  initial  velocity  of  0.01  mm/^ts  (10  m/s).  This 
condition  could  be  compared  to  a  particle  being  released  into  a  zero-temperature  bath  of  IV  —  1 
particles. 


1.3.5  Numerical  Approach 

The  original  code  was  written  by  Pfannes80  and  is  documented  in  appendix  A.  Modifications  were 
needed  to  support  studies  in  this  report.  The  most  significant  of  those  were  in  regard  to  the  decorated 
tapered  chain  (chapter  3)  and  is  documented  in  appendix  B.  In  these  numerical  studies,  the  velocity- 
Verlet  algorithm3  is  used  to  solve  the  differential  equations.  The  position  and  velocity  information 
are  updated  according  to: 


x(t  +  At.)  =  x(t)  +  v(t)At  H — —  a(t) 
v(t  +  At)  =  v(t)  +  ^-\a(t)  +  a{t  +  At)], 
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where  At  represents  the  timestep.  Accelerations  are  calculated  from  Newton’s  law  where  the  force 
is  evaluated  based  on  the  amount  of  overlap  determined  by  x(i),x(i  +  1).  Most  simulations  were 
evaluated  over  a  system  time  of  1  ms  where  the  timestep  was  set  to  10  picoseconds  with  108  steps  in 
the  integration  loop.  In  some  cases,  the  simulation  time  was  extended  to  10  ms.  For  each  simulation, 
a  separate  data  file  is  written  for  each  particle  which  contains  the  position,  velocity,  acceleration, 
kinetic  energy,  and  other  variables  per  timestep.  In  addition,  a  global  file  containing  the  potential, 
kinetic,  and  total  energy  of  the  whole  system  per  time  step  is  also  written  out. 

To  handle  the  thousands  of  simulations  that  needed  to  be  run,  an  automation  script  using  PERL 
was  written.  Appendix  C  lists  the  code  that  iterates  through  nested  loops  of  the  relevant  chain 
parameters.  In  each  cycle,  it  creates  new  subdirectories  based  on  the  current  value  in  the  loops, 
copies  a  template  containing  the  C++  source  listed  in  appendices  A  or  B  to  that  location,  replaces 
the  parameters  with  the  current  values  in  the  loop,  compiles  the  code,  and  then  runs  it.  This 
is  iterated  2200  times  for  simulations  comprising  chapter  2  and  about  1000  times  for  simulations 
comprising  chapter  3. 

1.4  Reduced  Problem 

In  order  to  gain  an  understanding  of  the  more  complicated  motions  for  systems  where  3  <  N  <  20, 
it  is  useful  to  look  at  single  and  binary  systems  confined  between  fixed,  but  compressible  walls.  It  is 
also  pedagogical  to  observe  the  changes  in  phase  space  when  one  moves  from  a  harmonic-like  (n  =  2) 
to  anharmonic  (n  yf  2)  dependence  in  the  overlap  potential. 

1.4.1  Single  Particle  System 

To  highlight  the  analytical  difficulties  in  dealing  with  equations  such  as  that  in  (1.6),  consider  the 
simplified,  but  general  problem, 

x  +  a:(n_1)  =  0,  (1.7) 

with  initial  conditions,  ;r(0)  =  0,x(0)  =  Vo,  where  n  >  2  and  values  may  be  fractions.  Surprisingly, 
little  work  has  been  done  in  dealing  with  solutions  to  second  order  differential  equations  where 
fractional  powers  need  to  be  considered.  There  have  been  several  efforts  by  Gottlieb88  and  the 
poineering  work  of  Mickens  6 '  who  have  used  the  method  of  harmonic  balance  which  uses  a  truncated 
Fourier  series  to  approximate  the  solution.  A  great  benefit  to  this  method  is  that  is  doesn’t  require 
perturbation  of  linear  terms,  so  it  can  —  in  principle  —  work  with  strongly  nonlinear  equations. 
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Figure  1.2:  Phase  space  diagrams  for  equation  (1.8). 


However,  efforts  have  so  far  have  failed  to  produce  a  closed  form  solution  to  x(t )  in  equation  (1.7). 
One  obtains  the  first  integration  of  equation  (1.7)  as, 


x  = 


(1.8) 


Without  explicitly  solving  the  EOM,  one  can  certainly  plot  the  phase  space  ( x ,  x)  for,  say,  n  = 
2.0,  2.5.  This  is  illustrated  in  figure  1.2  for  a  variety  of  vq  (arbitrary  units).  Starting  at  v(x  =  0) 
in  each  plot,  compression  increases  as  the  particle  slows  down.  For  r>o  =  1,  the  Hertz  potential  is 
soft  allowing  greater  compression  than  that  for  n  =  2.  This  was  also  reported  in  figure  1.1.  As  n 
increases  the  nonlinear  effects  become  apparent  as  the  Hertz  potential  becomes  increasingly  more 
repulsive  with  compression.  One  should  expect  therefore,  a  phase  space  similar  to  the  above  when 
looking  at  the  numeric  solution. 

To  calculate  the  period  of  oscillation,  T,  note  that  in  a  quarter-cycle  the  maximum  compression 
xm  corresponds  to  v  =  0.  The  period  can  then  be  written  as, 


T  = 


dx 


(1.9) 


For  n  =  2,  one  obtains  the  trivial  solution,  T  =  2ir.  It  is  clear  that  for  N  such  equations  coupled  by 
particle  overlaps,  one  needs  to  resort  to  numerical  methods  to  evaluate  Xi(t). 
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Figure  1.3  highlights  the  numeric  results  for  the  harmonic-like,  n  =  2  and  Hertz,  n  =  2.5,  cases. 
Both  are  oscillators  and  what  stands  out  immediately  is  the  “softness”  of  the  n  =  2.5  potential. 
This  is  visible  both  in  panels  (a) — where  in  the  anharmonic  case  the  period  is  longer — and  (b) 
where  the  closed  trajectories  in  the  anharmonic  case  clearly  reach  larger  displacements  in  x.  Panel 
(c)  demonstrates  the  relative  strengths  of  F  as  a  function  of  the  displacement.  It  also  appears 
that  the  initial  velocity,  Vi  =  lOm/s  does  not  cause  a  compression  significant  enough  for  the  Hertz 
potential  to  become  stronger  than  a  quadratic.  As  such  the  sphere  always  behaves  as  a  soft  particle 
whose  softness  increases  with  n.  Note  that  these  results  are  consistent  with  expectations  from  figure 

1.2  where  vq  =  1. 

One  can  perform  a  quick  study  on  a  single  particle  in  an  overlap  potential  well,  Sn.  In  this 
particular  case,  a  soft  particle — barely  touching  the  boundaries — has  an  initial  velocity,  ty.  One 
can  ask  various  questions  such  as  how  does  the  period  of  oscillation  scale  with  n  and  More 
importantly,  is  it  possible  to  obtain  a  single  expression  that  evaluates  T  =  T(n,Vi)l  Figure  1.4 
shows  the  decay  in  period  for  increasing  Vj  as  well  as  the  relative  increase  with  n.  The  latter  should 
be  expected  since  the  “hardness”  of  the  potential  is  proportional  to  n.  Because  of  the  excellent  fit 
afforded  by  a  power  law,  it  is  suspected  that  a  single  formula  can  describe  the  solution  space.  As 
such  we  write, 

T(n,Vi)  =  A(n)v~x(-n\  (1.10) 

where  A(ri),x(n)  can  be  fit  by  some  yet-to-be-determined  function.  It  turns  out  that  Vi  decays 
as  simply,  (n  —  2)  / n  while  the  values  comprising  A  decay  exponentially  with  a  single  phase  (3  and 
plateau  at  7.  These  conclusions,  including  the  fitting  functions  and  associated  coefficients  are  plotted 
in  figure  1.5.  This  then  allows  us  to  write  the  period  of  a  particle  oscillating  in  an  overlap  potential 
of  exponent  n  with  initial  velocity  Vi  as, 

T/2  =  (ae~0n  +  ~t)v~{n~2)/n  (1.11) 

where  the  coefficients  are  identified  in  the  figures  and  for  n  =  5/2,  T  oc  v-1/5. 

1.4.2  Two-Particle  System 

Figure  1.6  illustrates  the  velocity  and  position  phase  spaces  of  the  binary  system  for  the  first  120  /is. 
Several  points  of  interest  are  identified  in  the  velocity  phase  space — along  with  their  corresponding 
location  in  position  space — and  represent  either  a  zero  crossing  or  occasional  extremum.  Figure 
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-  n=2.0  ::  Harmonic  (Hooke) 

n=2.5  ::  Anharmonic  (Hertz) 


Figure  1.3:  Single  particle  system  between  fixed  compressible  walls  under  harmonic  (n  =  2.0)  and 
anharmonic  (n  =  2.5)  potentials  (a)  normalized  kinetic  energy,  (b)  phase  space,  (c)  force  versus 
displacement. 
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Figure  1.4:  Half-periods  plotted  as  a  function  of  v,  for  various  values  of  n.  Each  superimposed  fit 
has  the  form  given  by  T/2  =  A(n)vi  x^n'>  (dashed).  For  each  n  the  values  of  A,x  are  plotted  in 
Figure  1.5. 


Figure  1.5:  Coefficients  x,A  from  Figure  1.4  plotted  vs.  n  and  fitted  with  appropriate  functions  in 
support  of  generalized  expression  for  T/2. 
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1.7  depicts  these  points  in  a  cartoon  which  qualitatively  reproduces  the  velocities,  positions,  and 
compression  or  overlap  of  the  the  head  and  tail  particle.  Note  that  negative  motion  implies  movement 
towards  the  left  of  figure  1.7. 

It  is  clear  that  the  dynamics  have  become  much  more  complex  just  by  adding  an  additional 
particle.  As  one  moves  from  position  A  to  B,  the  head  particle  gives  up  its  energy  to  the  tail, 
becoming  maximally  compressed  at  position  C  in  panel  (b)  —  denoted  by  a  zero  crossing  in  panel 
(a).  Nearly  vertical  and  horizontal  lines  in  panel  (a)  denote  one  particle  accelerating  rapidly  while 
the  other  travels  at  a  small,  and  nearly  constant,  velocity.  The  only  way  this  can  happen  is  if  the 
particles  lose  contact.  That  is  the  overlap  vanishes,  <5  <  0.  These  instances  are  visible  in  figure  1.8 
which  plots  the  overlap  as  a  function  of  time  during  the  interval  seen  in  Figures  1.6,  and  1.7.  The 
ordinate  is  normalized  by  10  mm  which  is  the  length  of  the  combined  particles’  radii.  Loss  of  contact 
is  clearly  seen  in  figure  1.7  for  cases  G,J,  and  the  start  of  F.  The  consequence  is  momentary  faster 
particle  velocities  as  there  is  no  opposition  to  the  motion. 
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Figure  1.6:  Velocity  and  position  phase  diagrams  for  the  two  particle  system.  Dynamics  are  shown 
for  the  first  120  /is  with  corresponding  points  of  interest  identified.  Note  that  maximum  velocities 
do  not  necessarily  correspond  to  crossing  the  t  =  0  positions. 
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Figure  1.7:  Cartoon  depicting  the  evolving  two  particle  system.  Refer  to  figure  1.6  for  locations  in 
phase  space. 
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Figure  1.8:  Overlap,  <5,  as  a  function  of  time  for  the  particle-particle  interface  in  the  N  =  2  system. 
S  is  represented  as  a  percentage  of  the  two  spheres  combined  radii  since  it  is  normalized  by  10  mm 
(2  Ri). 
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Chapter  2 


The  Simple  Tapered  Chain  (STC) 


2.1  Introduction 

This  chapter  presents  an  in-depth  and  greatly  extended  study  of  TCs  originally  carried  out  by 
Pfannes80.  The  results  and  sections  of  this  chapter  were  reported  by  Doney  and  Sen25,26  and 
presented  at  the  2005  American  Physical  Society’s  topical  meeting  of  shock  compression  in  condensed 
matter  in  Baltimore,  MD.  Specifically  this  work  focuses  on  the  ability  to  spread  impulses  out  in  time 
and  space  rather  than  analyzing  solitary  wave  propagation  which  is  typically  the  dominating  topic. 
Section  2.2  addresses  the  hard-sphere  approximation,  with  and  without  a  term  accounting  for  energy 
loss  at  each  collision.  Section  2.3  looks  at  the  numeric  solution  in  terms  of  temporal  behavior  as 
well  as  normalized  kinetic  energy  (KE)  and  force  parameter  spaces.  That  section  also  includes  a 
discussion  of  how  these  systems  partition  energy.  This  is  followed  by  a  generalization  of  the  systems’ 
characteristic  fluctuations  and  the  approach  to  a  state  referred  to  as  quasi-equilibrium.  The  chapter 
concludes  with  a  mathematical  fit  to  describe  the  normalized  energy  mitigation  as  a  function  of  q,u>. 

2.2  Hard-Sphere  Approximation 

Hard-sphere  approximations  differ  from  the  numerically  simulated  systems  in  two  major  ways.  The 
first  is  that,  for  the  former,  the  chain  is  not  bounded  by  fixed  rigid  walls.  As  a  result,  energy  will  not 
continue  to  be  transmitted  up  and  down  the  chain.  The  second  is  that  the  potential  becomes  infinite 
( n  — >  oo )  and  as  a  consequence,  the  energy  packet  is  only  1  grain  in  width.  The  system  therefore 
propagates  energy  as  independent  collisions  and  is  congruent  to  the  model  proposed  by  Wu121. 
By  generating  an  iterative  form  of  the  conservation  equations,  one  can  arrive  at  an  expression  for 
the  normalized  kinetic  energy  (KE),  Ek  =  KEout/ KEin.  This  ratio  will  be  the  primary  variable 
determining  the  absorptive  quality  of  TCs. 

The  simple  tapered  chain  (STC)  is  displayed  in  figure  2.1.  To  generate  an  initial  disturbance,  an 
input  velocity,  ry,  is  applied  to  the  rightmost  and  largest  grain  with  radius  r,.  It  propagates  to  the 
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Figure  2.1:  Simple  tapered  chain:  N  =  10,  qs  =  8%,  L  =  70.7  mm,  r*  =  5  mm.  The  rightmost 
particle  is  grain  i,  its  nearest  neighbor  is,  i  +  1,  etc. 

left,  encountering  an  (initially)  stationary  grain  of  radius,  .  The  radius  of  the  (*  +  1)  particle 
may  be  reduced  to  (1  —  q)ri.  This  tapering  q  will  be  constant  along  the  entire  length  of  the  chain. 
During  the  transmission  of  the  impulse  along  the  chain,  there  may  be  energy  losses  and  we  consider 
two  cases  described  in  the  following  subsections. 


2.2.1  Lossless  STC  Hard-Sphere  Approximation 

Ignoring  any  energy  loss  during  a  collision,  we  perform  a  STC  hard-sphere  approximation.  Masses 
and  radii  are  expressed  as 


n+ 1  =n  -  rtqs  =  (1  -  qs)ri  =  eri 

mi  =  pVi  =  -7T  rfp  =  (2.1) 

mi+i  =  Wi+i  =  ?7£3rf ,  (2.2) 

where  e  =  1  —  qs.  Evaluating  the  conservation  of  momentum,  with  a  single  prime  denoting  post 
collision  values  and  the  initial  condition  that  the  (i  +  1)  particle  is  stationary  before  a  collision 
{vi+ 1  =  0),  all  rj  cancel  and  we  obtain 


rriiVi  +  mi+ivi+i 


Vi 


miv'i  +  mi+  n4+1 

rivi  +  1 

v'i  +  e3u(+1, 


(2.3) 
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where  equations  (2.1,  2.2)  have  been  used.  Following  the  same  procedure  for  the  conservation  of 
energy  while  ignoring  the  factor  of  one-half  yields, 


2  /2 
Vi  =  V. 


3„,/2 


=  v'f  +  e°v; 


i+i- 


(2.4) 


Letting  A  =  e3v'i+1  we  can  rewrite  equation  (2.3)  in  terms  of  v\  and  substitute  the  resulting  expression 
into  equation  (2.4), 

vi  —  {vi  —  A)2  +  Av'i+l 

=  v'f  ‘lAvi  +  A 2  T  Avi+i 

2  Vi  =  A  +  v'i+1 


L'i+i 

Vi 


1  +  e: 


3  ’ 


(2.5) 


Note  that  for  one  collision 

KEout  =  I<E'i+1  =  ml+ 1  (  v'i+l  \  2  =  g3/  v'i+1  \  2 
KEin  KEi  rrii  \  Vi  )  \  i>i  J 


KE'+1  4e3 

(1  +  e3)2  ' 


(2.6) 


For  N  particles  there  will  be  N  —  1  collisions,  each  of  which  have  the  ratio  in  equation  (2.6). 
Therefore,  the  normalized  KE,  Ek,  for  the  lossless  STC  hard-sphere  approximation  is  given  as 


Ek 


4[1  ~  # 

(l  +  [l-g]3)2/ 


(2.7) 


2.2.2  Lossy  STC  Hard-Sphere  Approximation 

The  same  approximation  can  be  performed  with  some  amount  of  energy  loss,  El,  per  collision. 
Consequently,  the  momentum  equation  (2.3)  is  unchanged,  but  equation  (2.4)  becomes 


=  v?  +  e3v'ii1  +  EL. 


(2.8) 
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One  then  obtains  a  more  complicated  expression  replacing  equation  (2.5): 


ui+ 1 
Vi 


2- 


Ejl 


1  +  e3 


(2.9) 


We  can  make  the  substitution,  El  oc  Viv'i+1  or  El  =  ELViv[+1,  where  El  is  the  constant  of 
proportionality.  This  adjustment  yields 

v'i+i  =  2e3  -  El 

Vi  e3(l  +  e3)  ’ 


The  corresponding  result  for  the  normalized  KE  for  N  particles  is 


2(1  —  q)3  —  El\ 2  \ 

(l-g)3[l  +  (l-g)3]2J 


(2.10) 


In  the  limit  El  =  0,  equation  (2.10)  reduces  to  the  lossless  case,  equation  (2.7),  as  one  would  expect. 
Note  that  results  are  independent  of  initial  velocity  and  size  of  the  grains. 


2.2.3  KE  Parameter  Space  for  STC  Hard-Spheres 

Figure  2.2  highlights  the  behavior  of  equation  (2.10)  for  0  <  q  <  0.1,  3  <  N  <  20  and  selected  El- 
Results  indicate  a  greatly  decreasing  output  of  energy  for  modest  values  of  El-  The  decay  of  N  also 
takes  on  an  exponential  decay.  Behavior  of  the  tapering,  q ,  resembles  sigmoidal  or  gaussian  decay 
and  approaches  linear  behavior  with  decreasing  N.  This  can  be  seen  in  equations  (2.7)  and  (2.10) 
and  is  expected  since  for  N  =  1  there  are  no  collisions  and  q  becomes  meaningless  (it  is  defined  for 
a  particle  and  its  neighbor) . 

Interestingly,  if  the  initial  velocity  is  supplied  to  the  smaller  end  of  the  TC,  one  also  observes 
shock  absorption  similar  to  the  system  in  figure  1,  albeit  with  less  efficiency.  In  this  case,  subsequent 
particles  are  growing  in  size,  i.e. ,  r,+ \  =  (1  +  g)r*.  Equations  (2.7)  and  (2.10)  are  modified  accord¬ 
ingly,  and  the  result  for  a  lossless  system  is  illustrated  in  figure  2.3.  Apparently,  both  configurations 
mitigate  a  propagating  pulse.  Note  that  the  range  of  axes  in  the  figures  are  different. 


2.3  Numeric  Solution 

Simulations  were  performed  for  3  <  N  <  20  and  0  <  q,  u>  <  0.1. 
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Figure  2.2:  Normalized  Ek{N,  q,  El)  parameter  space  for  the  STC  hard-sphere  approximation.  N 
varies  from  3  to  20  and  q  from  0%  to  10%. 
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Figure  2.3:  Normalized  Ex{N,q)  parameter  space  for  the  STC  hard-sphere  approximation  where 
initial  velocity  is  supplied  to  the  smaller  end  of  the  chain. 

2.3.1  Temporal  Behavior 

Figures  2.4,  2.5  sketch  the  temporal  behavior  of  the  smallest  grain  for  several  q  in  a  STC  where 
u>  =  0.05  and  N  =  15,  20,  respectively.  For  clarity,  these  plots  only  represent  the  first  half  of  the 
total  simulation  time.  Note  the  decreasing  scale  of  kinetic  energy  in  each  underlying  panel  in  both 
figures.  This  may  initially  appear  counterintuitive:  because  KE  scales  as  v2  and  m1  one  might 
expect  it  to  increase  with  larger  tapering  due  to  the  higher  velocities.  However,  m  oc  r3  ~  (1  —  q )3 
and  this  dominates.  Width  of  the  peaks  —  which  is  related  to  particle  velocities  —  are  functions  of  q 
and  N.  As  q  increases,  the  initial  peaks  are  shifted  to  earlier  times  so  signal  transmission  is  therefore 
faster.  Higher  velocities  also  imply  more  collisions  and  therefore  that  many  more  restitutive  phases. 

Note  that  in  both  cases,  a  single,  well-defined  pulse  has  been  turned  into  low  amplitude  noise 
by  increasing  the  tapering  in  the  system.  To  better  understand  the  dynamics,  we  can  observe  the 
behavior  of  several  juxtaposed  grains  near  the  boundary. 

Figure  2.6  highlights  the  absolute  positions,  velocities,  and  kinetic  energies  for  the  last  five  grains 
(particles  11-15)  in  a  chain.  In  particular,  let  us  examine  the  chains  represented  by  the  top  and 
bottom  subplots  of  figure  2.4:  monodisperse  and  10%  tapering,  respectively.  Particle  15  is  the  last 
grain  and  is  in  contact  with  the  boundary.  The  subplots  on  the  left  represent  the  monodisperse 
chain  where  each  grain  has  r  =  5  mm.  The  right  half  corresponds  to  a  chain  with  a  tapering  of 
10%  so  that  particles  11-15  have  radii  1.74,  1.57,  1.41,  1.27,  and  1.14  mm,  respectively.  Negative 
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Figure  2.4:  KE  as  a  function  of  time  for  the  smallest  particle  in  a  chain  with  N  =  15,  ui  =  0.5.  Each 
plot  is  evaluated  for  a  different  tapering  with  initial  kinetic  energy  as  0.0838  J. 
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Figure  2.5:  KE  as  a  function  of  time  for  the  smallest  particle  in  a  chain  with  N  =  20,  u>  =  0.5.  Each 
plot  is  evaluated  for  a  different  tapering  with  initial  kinetic  energy  as  0.0838  J. 
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Figure  2.6:  Absolute  positions,  velocities,  and  kinetic  energies  for  particles  11-15.  Particle  15  is  the 
last  grain  and  is  in  contact  with  the  boundary.  The  subplots  (a-c)  represent  the  monodisperse  chain 
where  each  grain  has  r  =  5  mm,  while  (d-f)  corresponds  to  a  chain  with  a  tapering  of  10%  so  that 
particles  11-15  have  radii  1.74,  1.57,  1.41,  1.27,  and  1.14  mm,  respectively.  These  subplots  illustrate 
when  the  grains  in  question  first  receive  the  incident  impulse.  Note  the  earlier  time  of  arrival  for 
q  =  0.1  since  velocities  are  higher  (up  to  a  factor  of  3  in  the  plots).  For  panels  (d-f)  the  dynamics 
have  been  extended  beyond  initial  incidence  and  reflection. 
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Figure  2.7:  Each  of  the  plots  depict  force  as  a  function  of  time  for  the  tail  and  head  particles  in  a 
chain  with  N  =  15,  u j  =  0.5.  The  value  of  the  tapering  increases  for  successive  plots  and  is  denoted 
by  the  label  q.  Negative  forces  imply  acceleration  towards  the  smaller  end  of  the  chain. 
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velocities  and  displacement  imply  motion  towards  the  end  of  the  chain. 

For  the  monodisperse  chain,  the  last  five  particles  first  receive  and  reflect  the  propagating  impulse 
between  40-120  /us.  It  also  appears  that  for  the  monodisperse  chain,  there  is  about  5  /us  of  compression 
for  a  given  grain  before  it  noticeably  responds  to  the  contact  force  and  accelerates  away.  Observe 
particle  11  acting  upon  particle  12  in  panel  (b)  for  example.  Note  that  after  another  5  /is  particle 
11  hits  its  peak  velocity  which  is  just  about  the  time  particle  13  starts  to  respond  to  the  impulse. 

When  looking  at  the  peak  KE  of  particle  13  occuring  at  62.5  /is,  energy  is  primarily  shared 
among  its  two  nearest  neighbors  such  that  most  of  the  kinetic  energy  of  the  system  is  shared  among 
3  adjacent  particles.  The  actual  number,  while  still  contested,  is  closer  to  5  or  7,  since  there  is  a 
small,  but  measurable  amount  of  energy  in  the  next  layer  of  nearest  neighbors.  This  is  true  for  any 
interior  particle  and  is  in  agreement  with  Manciu,  et.  al.59,60.  These  solitary  waves  (SW)  —  the 
“envelope”  of  the  five  grains  in  question  -  -  represent  localized  energy  that  maintain  their  width  yet 
decrease  in  amplitude  after  colliding  with  a  boundary.  The  decrease  of  the  tail  particle’s  velocity  is 
a  result  of  the  modified  potential  at  the  boundary. 

At  approximately  77  /us,  nearly  all  energy  is  briefly  stored  as  elastic  potential  energy.  Simulta¬ 
neously,  particles  14  and  15  are  maximally  displaced  from  their  equilibrium  position  while  particles 
11-13  are  still  slowly  compressing  towards  the  boundary.  A  close  examination  about  this  time  shows 
that  the  fourteenth  grain  is  the  first  particle  to  reverse  direction  due  to  the  restoring  force  of  the 
last  particle.  Acceleration  plots  indicate  more  clearly  that  this  happens  while  the  tail  particle  is  still 
being  compressed  with  the  boundary. 

The  big  picture  of  energy  propagation  then  consists  of  several  regimes.  In  the  microscopic  sense 
one  could  describe  the  elastic  waves  emanating  from  the  contact  region  and  reflecting  within  any 
particular  grain.  Mesoscopic  features  have  been  the  current  focus  of  these  plots  and  observations  — 
that  is,  bulk  motion  of  a  grain.  And  the  macroscopic  picture  would  be  the  solitary  (longitudinal) 
wave  formed  by  several  grains  and  discussed  in  previous  work. 

For  the  highly  tapered  chain  (q  =  0.1)  on  the  right  hand  side  of  Figure  2.6,  the  last  five  particles 
first  receive  and  reflect  the  propagating  impulse  in  about  half  the  time  as  the  monodisperse  case.  As 
particles  11-14  are  accelerated  towards  the  end,  their  maximum  speed  increases  due  to  their  smaller 
masses.  Particle  11  accelerates  at  about  25  /us  and  compresses  particle  12  for  nearly  2  /its  before  it 
starts  to  move.  Shortly  thereafter,  particle  11  slows  down  at  about  28  /us  due  to  the  contact  force 
with  its  neighbor,  particle  12,  and  additional  contact  forces  with  particles  13-15.  This  is  also  true 
for  down  range  particles  12  and  13;  however,  the  accelerations  (not  pictured)  are  enhanced  as  one 
moves  closer  to  the  boundary. 
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What  stands  out  dramatically  in  panel  (f)  is  the  absence  of  complete  energy  transmission.  All 
particles  shown  have  a  measurable  amount  of  energy  when  the  tail  particle  rebounds  at  about  32.5  /is 
which  is  absent  in  the  monodisperse  chain.  The  displacement  plot  shows  all  grains  biased  towards 
the  tail  of  the  chain  from  the  equilibrium  position.  Each  minima  corresponds  to  a  maximum  in 
their  elastic  potential  energy.  Because  these  occur  at  slightly  skewed  times,  subsequent  collision- 
compression  oscillations  occur  which  are  much  more  visible  in  the  following  50  fxs  (not  pictured).  It 
is  difficult  to  draw  many  conclusions  from  such  a  complicated  landscape.  However,  it  is  expected 
that  no  solitary  wave  is  present  in  such  chains:  the  dispersion  due  to  tapering  is  too  great. 

When  one  records  the  force  or  energy  at  the  end  of  the  tapered  chain  with  a  sensor  of  some 
sort,  the  patterns  recorded  vary  for  N  and  q  as  these  plots  suggest.  As  the  tapering  increases,  the 
pattern  changes  from  a  well-defined  and  periodic  pulse  to  noise.  This  trend  is  visible  for  increasing 
q  in  Figures  2.4,  2.5  and  represents  the  thermalization  that  we  are  looking  to  exploit  in  various 
applications. 

2.3.2  KE  and  F  Parameter  Space  Behavior 

The  effectiveness  of  TCs  can  be  measured  based  on  the  normalized  kinetic  energy  Ek  and  force 
surfaces  formed  by  varying  N,  q ,  and  u>,  thereby  creating  many  tapered  chain  configurations.  Specif¬ 
ically,  we  form  the  ratio,  Ek  =  KEout/KEin,  where  KEout  is  the  first  peak  felt  by  the  last  grain 
and  KEin  is  unchanged.  Analogously,  Fout  is  the  first  minimum  felt  by  the  last  grain  since  the 
direction  of  a  negative  force  is  into  the  wall  or  force  sensor.  The  algorithm  to  pick  out  the  first 
turning  point  for  each  of  these  is  straightforward.  Since  KE(t )  is  a  column  vector,  one  can  iterate 
through  each  element  until  the  first  occurence  of  (i  +  1)  <  i  is  true.  In  that  case,  i  represents  the 
peak.  For  force  we  compare  against  (i  +  1)  >  *.  Actually,  normalization  is  a  critical  component  in 
determining  the  effectiveness  of  TCs  and  has  been  given  a  more  thorough  discussion  in  Appendix  E. 

Figure  2.8  highlights  the  numerical  results  for  KE(N,  q ,  constant  u>).  The  gaussian  and  exponen¬ 
tial  dependence  on  q  and  N,  respectively  stand  out  immediately.  These  KE  surfaces  represent  STC 
chains  that  thermalize  more  than  half  the  incident  energy  introduced  into  the  system.  For  example, 
the  least  effective  shock  mitigating  geometry  that  we’ve  simulated  here  reduces  the  output  kinetic 
energy  to  approximately  40%  of  the  input  which  is  shown  in  panel  (a)  where  N  =  3,g  =  0.  The 
“missing”  energy  is  distributed  among  the  other  grains  as  kinetic  and  potential.  Note  that  results 
are  independent  of  initial  grain  size.  For  monodisperse  chains  (q  =  0),  there  is  asymptotic  behavior 
as  N  increases.  That  is,  when  u>  =  0  for  monodisperse  chains,  KE  becomes  a  constant.  This  result 
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is  consistent  with  observations  by  Sen,  et.  al. 96  that  approximately  15  grains  are  required  for  the 
solitary  wave  to  become  fully  established.  This  then  propagates  indefinitely  (in  our  ideal  setup)  so 
that  additional  particles  for  to  =  q  =  0  have  little  to  no  effect.  As  restitution  is  adjusted  to  higher 
values,  the  decompressive  force  is  reduced  to  (1  —  lo)%  of  the  compressive  force.  Consequently,  even 
the  monodisperse  chains  become  lossy  for  increasing  N. 

For  small  chains  there  are  few  collisions,  so  restitutive  losses  do  not  play  a  very  big  role.  The 
same  is  true  for  the  amount  of  tapering.  If  one  holds  q  fixed  at  a  small  value  and  increases  N, 
the  chain  is  almost  monodisperse.  The  mismatches  are  insignificant  enough  that  localized  energy 
can  propagate  and  only  over  large  N  does  it  break  up.  This  translational  symmetry  breaking  was 
also  reported  by  and  is  in  agreement  with  Nakagawa,  et.  al'2.  Coupled  with  restitutive  losses, 
the  longest  and  most  tapered  chains  quickly  lose  energy  through  the  massive  number  of  collisions 
occuring.  As  a  shock  mitigating  technology  however,  one  wants  to  maximize  the  thermalization  even 
for  small  chains  since  there  may  not  be  the  luxury  of  large  volumes  in  certain  applications  (hence 
the  improved  chain  design  presented  in  chapter  3. 

A  similar  surface  can  be  plotted  representing  the  normalized  force  for  various  tapered  chains. 
This  is  visible  in  Figure  2.9  and  is  similar  to  the  corresponding  KE  plots.  Again,  the  functional 
form  of  N  appears  as  a  one  or  two  phase  exponential,  but  the  gaussian  nature  of  q  is  less  obvious 
and  can  be  approximated  linearly  throughout  the  parameter  space.  In  general,  the  magnitude  of 
TV  appears  to  be  twice  that  of  Ek- 

The  KE  numerical  results  can  be  compared  to  the  hard  sphere  approximation  by  plotting  their 
difference  (Figure  2.10).  Comparisons  that  include  losses  are  not  evaluated.  One  difference  between 
the  two  parameter  spaces  that  stands  out  occurs  for  small  N.  As  one  traces  out  a  path  along 
constant  q  =  0  and  increasing  N,  the  numerical  results  (Fig.  2.8a)  show  a  quick  decrease  in  Ek 
which  then  approaches  an  asymptote.  We  suspect  that  this  is  due  to  the  width  of  input  pulse  being 
longer  than  the  distance  to  the  boundary.  This  does  not  occur  in  the  hard-sphere  approximation 
since  independent  collisions  preclude  that  possibility. 

2.3.3  Energy  Partition 

The  partitioning  of  energy  in  a  multi-bodied,  dynamical  system  is  a  helpful  tool  in  understanding 
the  behavior  of  discrete  systems.  For  instance,  it  provides  a  measure  of  how  much  of  the  system  is 
in  motion  versus  its  response  to  the  potential.  Energy  information  about  the  whole  TC  system  can 
be  evaluated  by  summing  the  kinetic  and  potential  contributions  of  each  bead.  In  all  cases,  we  have 
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Figure  2.8:  Numerical  solution  of  Ek(N,  q)  parameter  space  for  constant  u> 
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Figure  2.9:  Numerical  solution  of  F(N,  q)  parameter  space  for  constant  u> 
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Figure  2.10:  Difference  plot  of  the  lossless  KE  parameter  spaces  for  the  hard-sphere  approximation 
(Figure  2.2a)  and  numerical  solution  (Figure  2.8a).  Note  that  the  azimuthal  view  has  been  rotated 
180  degrees  for  clarity. 
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Figure  2.11:  Energy  partitioning  for  STCs  where  N  =  {3,  8, 14,  20}  represents  the  number  of  spheres 
and  qs  =  {0.0,0.05,0.1}  is  the  tapering.  Noisy  plots  resembling  panel  l  indicate  shock  absorbing 
systems  while  those  similar  to  panel  j  transmits  an  impulse  as  a  solitary  wave  —  essentially  without 
loss.  Panels  a-c  represent  efficient  energy  conversion  systems. 


seen  that  energy  sharing  among  particles  is  rapid  and  remains  nonlocalized.  Here  attention  is  paid  to 
the  energetic  response  of  the  system  with  some  consideration  to  individual  particle  velocities.  Note 
that  q  is  now  being  written  as  qs  to  distinguish  it  from  results  that  will  be  presented  in  Chapter  3. 
Therefore  q  =  qs  represents  the  tapering  in  a  STC. 

Energy  partitioning  in  the  STC  is  displayed  as  twelve  subplots  in  figure  2.11  where  each  plot 
element  represents  a  chain  with  N  =  {3,  8, 14,  20}  and  qs  =  {0,  0.05,  0.1}.  For  clarity,  only  the  first 
fifth  of  simulation  time  is  shown.  Many  interesting  features,  visible  in  the  plots,  can  be  discussed 
qualitatively  where  it  appears  that  simple  chains  can  be  loosely  categorized  into  three  energy  regimes: 
solitary  wave  systems  (small  qs,  large  N),  shock  absorption  systems  (large  qs ,  large  N),  and  strongly 
oscillating  systems  (small  N) .  For  smaller  systems,  it  is  still  unclear  how  exactly  the  system  behaves 
when  the  width  of  the  solitary  wave,  wsi  is  smaller  than  the  size  of  the  system,  N. 

The  simplest  case  to  consider  is  that  of  panel  (d)  as  it  has  been  presented  before80  and  represents 
a  monodisperse  chain  of  20  grains.  It  is  clear  in  this  instance  that  energy  is  partitioned  into  55-57% 
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Figure  2.12:  Normalized  KE  for  N  =  3,  qs  =  0.05.  Ei  represents  the  head  of  the  chain  and  input;  E2 
is  the  central  bead;  and  E3  is  the  tail.  Note  the  periodicity  of  the  system  resulting  in  a  recurrence 
time  of  just  under  80  /rs. 

kinetic  and  45-43%  potential.  An  inspection  of  grain  speeds  reveal  that  the  first  several  particles 
retain  some  residual  velocity  before  a  SW  can  form  —  which  also  holds  true  for  every  plot  in  the 
figure. 

Holding  N  =  14  fixed  and  increasing  qs  results  in  faster  particles  as  one  moves  down  the  chain 
as  well  as  wave  broadening.  This  is  visible  in  panels  (g,h,k,l)  as  an  increasing  ramp  whose  slope  is 
steeper  with  increasing  qs  and  smoothening  of  the  sinusoidal  modulations,  respectively.  Compare  this 
with  the  wave  broadening  among  multiple  grains  with  increasing  qs  seen  in  figure  2.13.  Interaction 
with  the  boundary  follows  and  in  some  cases,  the  energy  envelope  increases  because  most  particles 
have  already  reversed  direction.  Close  inspection  reveals  that  trailing  particles  can  catch  up  and 
kick  leading  particles  into  higher  energy  states.  In  panels  ( k,l )  rapid  oscillations,  or  thermalization, 
is  visible.  This  is  an  indication  of  increasing  “randomness”  111,119  of  motion  and  spreading  of  the 
energy  in  time  and  space  —  a  prerequisite  for  impulse  decimation. 

As  one  moves  to  shorter  chains,  where  ws  >  N  like  those  in  panels  (a,e,i),  nearly  complete  con¬ 
version  of  energy  frequently  takes  place.  Wave  reflection  therefore  has  begun  before  full  transmission 
of  the  incident  pulse  would  normally  reach  the  boundary  of  a  longer  chain.  Increasing  qs  for  N  =  3 
appears  to  enhance  the  energy  conversion  efficiency.  In  fact,  panel  (6)  is  particularly  alluring  given 
its  quasi-periodicity.  Amplitude  of  the  intermediate  weak  peaks  vary  based  on  qs  and  represent 
energy  transfer  through  the  center  bead.  This  behavior  throughout  the  simulation  suggests  the 
possibility  of  nonlinear  modes115.  Figure  2.12  provides  a  glimpse  into  future  work  on  the  subject 
where  the  normalized  KE  of  each  particle  is  superimposed  as  a  function  of  time.  The  recurrence 
time  —  or  period  of  the  system,  r  —  is  just  under  80  fj,s7  and  decreases  slightly  with  every  cycle.  It’s 
not  clear  whether  there  exists  a  critical  value  qc  that  admits  a  constant  recurrence  time.  Increased 
values  of  qs  tend  to  attenuate  E\.  Note  also  the  lack  of  equipartition  for  this  system. 
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Instantaneous  KE  per  grain  -  t=52us 


grain  number 

Figure  2.13:  Instantaneous  kinetic  energy  per  grain  for  N  =  15,  t  =  52/rs  and  selected  tapering. 

Additionally,  one  can  take  a  mesoscopic  view  and  investigate  how  STC  systems  break  down  the 
energy  per  grain.  Observe  figure  2.13  where  the  instantaneous  kinetic  energy  is  plotted  per  grain 
for  N  =  15,  t  =  52/rs  and  a  selection  of  tapering  values.  In  these  plots,  the  incident  impulse  moves 
from  right  to  left.  Panel  (a)  represents  the  monodisperse  chain  and  the  localization  of  energy  is 
apparent  as  the  SW  is  constructed.  Why  does  it  take  so  many  particles  before  a  SW  is  constructed? 
Note  from  Figures  2.11(b-d)  that  it  takes  a  finite  amount  of  time  —  or  equivalently,  distance  —  for 
the  system  to  partition  energy  according  to  the  virial  theorem  (see  section  2.3.4)  which  for  n  =  2.5 
implies  <  I\E  >=  0.55.  This  means  that  even  for  inertially  matched  spheres,  some  residual  energy 
is  left  at  the  head  of  the  chain.  Panels  (b,c)  again  illustrates  the  effect  of  tapering  which  spreads 
the  energy  out  among  grains.  For  the  latter,  the  distribution  before  interacting  with  a  boundary 
appears  geometric.  It  is  clear  that  a  SW  cannot  exist  in  such  chains  because  of  the  large  dispersion. 
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An  even  more  illuminating  series  of  images  can  be  generated  which  shows  the  energy  landscape 
for  all  particles  in  the  chain  at  all  points  in  time.  This  is  illustrated  in  figure  2.14.  This  series  of  plots 
highlights  the  propagation  of  a  pulse  in  time  (horizontal  axis)  through  the  chain  of  grains  (numbered 
vertically)  for  several  values  of  tapering,  qs  =  {0,0.02,0.04,0.06,0.08,0.1}.  In  each  case,  N  =  20. 
Kinetic  energy  is  displayed  and  has  been  magnified  several  times  in  order  to  pick  out  the  smaller 
amplitude  backscattering  and  wave  phenomena.  For  the  monodisperse  chain,  one  can  clearly  see  the 
SW  propagate,  with  negligible  loss  and  backscattering  for  these  early  times.  It  is  easy  to  see  that 
the  SW  velocity  is  approximately,  1  mm / /is,  this  is  in  contrast  to  the  delta  function  impulse  in  grain 
number  20  of  10  mm//us  —  recall  that  some  residual  energy  is  left  over  in  the  first  several  grains. 
This  velocity  decreases  ever  so  slightly  per  cycle  as  the  SW  is  destroyed  and  recreated  imperfectly 99 
at  the  boundaries.  As  q  increases,  trapped  energy  is  seen  by  the  light  horizontal  bars  in  panels 
(e,f)  preceding  100  [is.  What’s  also  very  interesting  to  note  here  is  the  increasing  magnitude  of 
backscattering  during  the  return  phase  to  particle  20.  With  increasing  q ,  energy  transfer  is  hindered 
during  the  momentum  transfer.  Starting  with  panel  (d)  at  100  /us,  many  subsequent  waves  of  lower 
velocity  towards  the  head  of  the  chain  are  visible.  At  500  /xs,  organized  waves  emerge  that  only 
propagate  through  the  smaller  half  of  the  chain  (top  half  of  panels) .  This  effect  also  occurs  for  larger 
q  in  panels  (e,f)  but  more  quickly  —  after  the  first  interaction  with  the  boundary.  In  these  last 
panels,  the  emergence  of  organized  additional  wave  motion  was  unexpected. 

2.3.4  The  Approach  to  Quasi-equilibrium 

What  is  the  equilibrium  state  for  a  nonlinear  system?  For  that  matter,  what  is  it  for  a  strongly 
nonlinear  system  with  no  restoring  terms  such  as  the  Hertz  potential?  These  questions  have  been 
addressed  by  previous  researchers  with  emphasis  on  FPU  systems  as  discussed  in  section  1.2.2.  Here 
we  outline  an  approach  that  leads  to  an  interesting  state  referred  to  as  quasi-equilibrium  (qeq) 
where  the  system  exhibits  some  traits  of  being  in  equilbrium.  It  is  found  that  (1)  the  eventual 
state  of  the  system  appears  to  be  initial  condition  dependent  —  this  is  opposite  the  expectation 
of  an  equilibrium  state;  (2)  the  distribution  of  velocities  of  the  system  is  roughly  a  gaussian  -  so 
a  maxwellian  distribution  of  velocities  is  found  —  this  is  expected  in  an  equlibrium  phase;  (3)  TC 
systems  appear  to  have  sustained  fluctuations  that  deviate  significantly  from  equilibrium. 
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Velocity  Distributions 


For  a  dilute  gas  in  equilibrium,  constituent  particles  exhibit  a  maxwellian  (or  gaussian)  velocity 
distribution84.  One  can  therefore  measure  whether  a  similar  system  has  reached  equilibrium  by 
binning  the  velocities  of  each  particle  and  observing  whether  the  profile  is  that  of  a  gaussian  —  as 
figures  2.15,  2.16  do  for  N  =  4,20,  respectively.  In  these  panels,  there  are  50  bins  ranging  from 
-0.01  to  0.01  mm / /js  where  the  vertical  axes  represents  the  count  per  bin.  Note  that  each  panel 
automatically  adjusts  the  number  of  counts.  In  both  cases,  particle  1  represents  the  end  of  the  chain 
and  particle  N  represents  the  head  where  the  initial  impulse  was  applied.  In  addition,  simulation 
time  was  extensive:  10  ms.  For  earlier  times,  the  profiles  differed  markedly  from  a  gaussian.  The 
plots  show  that  with  fewer  particles,  it  is  more  common  to  find  higher  velocities,  that  is  they  don’t 
decay  significantly.  For  large  N  however,  velocities  decay  to  the  point  of  a  super-abundance  near 
the  mean  value  —  note  the  large  increase  of  counts.  It  was  expected  that  boundary  particles  would 
stand-out  in  their  profiles;  however,  interior  and  boundary  particles  appear  indistinguishable  in  this 
representation.  Something  that  is  missing  from  these  plots  is  the  temporal  dependence  —  that  is, 
how  long  does  it  take  to  reach  such  profiles? 

Equipartition 

The  next  issue  that  can  be  addressed  is  whether  or  not  one  has  equipartition  of  energy.  The  virial 
theorem37,4'-103  describes  the  average  kinetic  and  potential  energy  in  a  system.  If  the  force  can  be 
derived  from  a  potential  of  degree,  n,  then  the  virial  theorem  becomes, 

-n  <  PE  >=<  KE  > 

2 

71 

<KE>=-& TT)  <211> 

where  the  normalized  total  energy,  E  =  1  =<  KE  >  +  <  PE  >,  has  been  used  to  eliminate 
<  PE  >.  This  is  the  expected  mean  value  of  KE  for  the  system.  Note  that  n  =  2  is  the  only  case 
where  the  energy  is  equally  divided  between  potential  and  kinetic.  Thus,  is  the  strictest  sense  of  the 
word,  equipartition  has  been  achieved.  One  can  compare  results  from  the  simulation  to  the  these 
expectations.  Figure  2.17  illustrates  the  average  KE  as  a  function  of  N  for  several  values  of  n.  It 
should  not  be  surprising  that  <  KE  >  is  independent  of  N  since  it  has  no  explicit  dependence  in 
the  Hamiltonian.  Figure  2.23  demonstrates  that  the  simulational  results  agree  with  what  one  would 
expect  from  equation  (2.11). 
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If  equation  (2.11)  is  equally  distributed  to  each  particle  in  the  system,  then  for  a  system  ap¬ 
proaching  equilibrium  one  might  expect 


<  KE  >n= 


n/(n  +  2) 
N 


(2.12) 


per  particle.  These  values  are  measured  in  figures  2.19,  2.20,  and  2.21  at  100,  500,  and  1000  /is, 
respectively.  Each  plot  represents  a  matrix  of  subplots  for  several  STCs  where  N  =  {3,  8, 14,  20}  and 
q  =  {0,  0.05,  0.1}.  The  horizontal  axis  in  each  panel  represents  the  specific  grain  and  the  vertical  is 
the  normalized  energy.  The  red  line  indicates  the  theoretical  expectation  given  in  equation  (2.12). 
As  the  sequence  is  observed  in  time  between  the  figures,  it  is  clear  that  short  chains  with  large  q 
do  not  equally  divide  energy  (recall  figure  2.12)  as  described  by  equation  (2.12).  Even  for  extended 
times,  longer  chains  do  not  always  reach  this  limit  either  —  particularly  for  grains  adjacent  to  the 
boundary  and  q  =  0.05  for  example. 


Fluctuations 

The  mean  value  of  KE  alone  doesn’t  say  much  about  the  system.  It  is  therefore  of  interest  to  measure 
the  fluctuations  about  this  mean  and  their  decay  with  increasing  N.  The  evaluation  of  fluctuations 
however,  is  usually  reserved  for  large  systems.  While  the  STCs  we  are  considering  only  have  up  to 
N  =  100  particles,  it  is  the  number  of  timesteps  in  the  simulation  that  is  actually  being  considered; 
and  that  quantity  is  large  with  50000  elements. 

The  fluctuations  represent  a  nearly  continuous  value  (dictated  by  the  size  of  the  time  step),  F(t)  = 
KE[t)—  <  KE  >.  To  keep  the  value  positive  definite  it  is  written  as,  F(t)  =  \J (KE(t)—  <  KE  >)2 
and  figure  2.22  illustrates  the  results  for  several  chains:  F(t,  n,N=  100,  q  =  0).  For  increasing  n  it 
appears  that  the  mean  value  of  the  fluctuations,  <  F  >,  is  growing.  In  addition,  the  time  to  reach 
the  boundary  is  increasing  which  is  denoted  by  the  peaks.  The  latter  observation  suggests  that  the 
spheres  become  increasingly  soft  with  n  which  is  counterintuitive  since  one  obtains  the  hard-sphere 
approximation  as  n  — »  oo.  What  is  happening?  The  answer  lies  in  Figure  1.3(c).  Recall  that  this 
figure  plots  the  force  versus  position,  or  compression,  for  a  single  grain  in  the  n  =  2,  2.5  potentials. 
Note  that  the  displacement  isn’t  large  enough  such  that  the  Hertz  force  overtakes  the  Hookean- 
like  force.  As  a  result,  for  increasing  n,  the  potentials  become  increasingly  soft  —  as  Figure  1.1 
indicates  —  for  small  displacements  (i.e.  those  less  than  1).  That  is  the  reason  why  for  increasing 
n,  interactions  appear  softer  even  though  theoretically,  one  is  approaching  the  hard-sphere  limit  — 
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recall  figure  1.2.  The  mean  value  of  the  fluctuations  can  be  calculated  by  the  following, 


<  F  >= 


NAt 

E 


i—At 


m 

NAt 


where  F)(f)  is  the  instantaneous  value  of  F,  At  =  0.2/rs  is  the  timestep,  and  NAt  =  50000  is  the 
number  of  timesteps. 

One  would  expect  the  fluctuations  to  approach  zero  as  N  tends  to  infinity47  since  there  are 
more  particles  to  distribute  the  energy  over  with  increasingly  smaller  changes  about  the  mean.  The 
current  data  however,  indicates  that  this  may  not  be  the  case  for  overlap  potentials.  It  is  found  that 
<  F  >  plateaus  at  a  value  that  is  a  function  of  n  —  see  for  example  figure  2.23(a).  Specifically,  for 
n  =  2  and  n  =  2.5  differences  between  data  is  small  when  one  looks  at  <  F  >  and  it  was  found 
that  in  both  cases,  <  F  >  is  about  the  same  and  decreasing  in  N.  Thus,  the  distinction  in  whether 
or  not  <  F  >— >  0  is  hard  to  make  when  n  =  2  and  n  =  2.5.  It  is  pretty  clear  however  that  when 
n  >  2.5,  <  F  >  appears  to  stabilize.  Thus,  these  systems  seem  to  have  sustained  fluctuations  that 
deviate  significantly  from  equilibrium.  The  reason  why  this  happens  is  because  the  solitary  waves 
that  carry  energy  span  several  particles  and  their  width  is  controlled  by  n.  Thus,  any  two  particles 
don’t  quite  carry  the  same  energy  on  average  at  any  given  time.  It  was  believed  then  that  when 
q  >  0,  <  F  >  would  approach  zero  since  one  is  forcing  thermalization  by  breaking  down  the  solitary 
waves.  This  is  not  quite  so  however  and  it  remains  a  topic  for  future  investigation. 

In  order  to  characterize  how  F  changes,  one  can  measure  <  F  >  as  a  function  of  n.  Therefore, 
figure  2.23(a)  plots  the  decay  of  <  F  >=<  F(n,  N )  >.  Each  set  of  data  appears  to  follow  a  function 
of  the  form, 

<  F  >=  A(n)e~k{n)N  +  B(n),  (2.13) 


where  Prism71  has  been  used  for  the  nonlinear  regression  analysis  in  determining  the  values  of 
A,  B ,  k  for  each  n.  Here,  A  represents  the  maximum  amplitude  above  the  saturation  B  while  k 
governs  the  phase  (or  half-life)  of  the  exponential  decay.  Note  based  on  the  trend  in  the  data,  this 
model  doesn’t  allow  <  F  >— >  0  in  the  limit  of  N  — »  oo.  Instead  it  saturates  at  B.  It  is  possible 
though  for  B  — >  0;  however  N  isn’t  sufficiently  large  for  that  to  potentially  occur.  For  n  =  2,2.5, 
<  F  >  appears  to  continue  to  drop,  although  our  maximum  value  of  N  would  need  to  be  increased 
to  prove  this  assertion.  For  the  range  in  data  however,  equation  (2.13)  fits  the  data  well. 

To  obtain  expressions  for  how  the  parameters  A1  B ,  k  vary  as  functions  of  n,  they  have  been 
plotted  in  figure  2.23(b-d).  In  panel  (b),  k  appears  scattered  and  independent  of  n.  For  simplicity, 
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we  therefore  will  use  its  mean  value,  <  k  >~  0.3.  Panels  (b,c)  however,  appear  to  follow  an 
exponential  decay  with  plateau  and  sigmoid,  respectively.  With  the  coefficients  selected  by  Prism 
as  a  best  fit  with  narrow  95%  confidence  intervals,  the  relations  are  given  as, 


A(ri) 

B{n) 


0.2052e_ol898ri  +  0.04246 

0.1563  +  0.07877 

—0.07877- 


1  +  exp  [(3.673  -n)/2.216]' 


(2.14) 


What’s  interesting  is  that  A(n  =  2)  and  B(n  =  2)  are  not  represented  well  by  the  fit  as  seen  in 
panels  (c,d).  This  is  odd  considering  that,  in  regard  to  B ,  the  data  points  along  the  plateau  in  panel 
(a)  for  n  =  2,  2.5  are  nearly  identical.  When  equations  (2.14)  and  <  k(n)  >=  0.3  are  combined  with 
equation  (2.13)  one  obtains  a  full  expression  for  the  decay  in  fluctuations  as  functions  of  n  and  N. 
The  results  are  plotted  in  Figure  2.24  where  values  for  n  —  2  have  been  excluded  since  they  fall 
outside  the  model’s  applicability  ( n  >  2.5,1  <  N  <  100).  The  data  is  well-fitted  by  the  model  and 
allows  one  to  interpolate  the  fluctuations  in  the  average  KE  of  a  system  governed  by  the  overlap 
potential,  Sn,  for  various  n  and  N.  Because  the  initial  condition  of  Vi  =  lOm/s  only  compresses 
grains  such  that  larger  n  leads  to  softer  grains,  it  would  be  useful  to  expand  (2.13)  to  take  into 
account  ry.  Note  that  the  values  of  F(t)  have  been  obtained  over  very  long  simulation  time,  104/zs, 
yet  the  fluctuations  about  <  KE  >  of  the  system  remain  measureable. 


Quasi-equilibrium 

For  an  isolated  system  to  be  in  equilibrium  it  implies  that  dynamic  properties  of  that  system  do  not 
change  with  time10' .  Take  for  example  <  KE  >,  it  was  seen  in  the  previous  section  that  while  some 
systems  may  reach  a  generalized  or  mixed  equipartition,  their  mean  values  of  <  I\  E  >  still  fluctuate 
in  time  quite  considerably.  This  was  quantified  in  equation  (2.13).  An  equilibrium  state  would 
require  the  fluctuations  to  vanish,  yet  we  observe  them  to  saturate  at  a  nonzero  value  even  over  long 
simulation  times.  When  one  additionally  takes  into  account  that  STC  chains  have  gaussian  velocity 
distributions,  it  seems  that  such  systems  display  some  properties  of  equilibrium,  but  not  all  of  them. 
This  curious  state  is  referred  to  as  quasi-equilibrium  (qeq)  and  is  a  very  recent  discovery 69,99,100 . 

One  way  to  characterize  this  state  is  to  measure  the  relaxation  time  to  qeq.  With  increasing 
N,  energy  spreads  out  and  consequently  so  do  the  maximum  velocities.  Rather  than  looking  at  the 
velocity  distributions  per  grain  therefore,  we  can  record  the  absolute  maximum  velocity,  \vrnax\  of 
any  grain  in  a  chain  per  timestep — or  over  some  interval  of  timesteps — and  plot  that  as  a  function  of 
time.  This  is  illustrated  in  figure  2.25  which  plots  \vmax\  as  a  function  of  time  for  N  =  16, 18,  20, 50. 
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Note  that  when  the  velocities  decay,  they  saturate  to  some  value.  The  relaxation  time  to  qeq,  Tqeq, 
therefore,  is  when  that  saturation  occurs  with  negligible  change  over  time.  Characterization  of  the 
time  is  difficult  however  because  it  is  calculated  from  a  fitting  function  based  on  data  which  has  a 
significant  amount  of  scatter.  In  a  first  attempt  to  measure,  Tqeq,  the  results  proved  inconclusive 
as  the  energy  relaxation  of  the  system  needs  to  be  better  characterized  as  to  reduce  the  amount  of 
scatter.  What  follows  is  a  description  of  the  technique  which  will  form  the  basis  of  future  work  on 
the  subject. 

The  decay  can  be  fitted  to  a  Gaussian  whose  peak  or  mean  is  centered  on  t  =  0  such  that  the 
fitting  function  is 

Y  =  .4  _  (2.15) 

where  A,B,a  are  variables  to  be  fit  (by  Prism);  however  in  each  attempt,  very  poor  residuals  were 
obtained  due  to  the  large  scatter  in  the  data.  The  fit  appeared  accurate  for  N  =  18,  20;  however  it 
was  less  than  certain  for  N  =  16,  50.  In  particular,  for  N  =  16,  there  is  only  a  small  difference  in 
t  =  0, 10000  values  —  the  chain  could  be  too  short  for  such  measurement. 

When  a  fit  is  established,  one  can  measure  rqeq  by  measuring  the  difference  in  subsequent  elements 
in  the  vmax  vector.  When  that  value  drops  below  a  certain  threshold,  the  remaining  differences  are 
negligible  and  the  system  has  reached  qeq.  This  is  expressed  as,  vmax{i  +  1)  —  Vmax(i)  <  vc  and 
based  on  the  plot  of  these  differences  as  a  function  of  time  (not  pictured),  one  can  chose  vc  =  1  ■  10”' 
mm//js.  Results  were  strongly  dependent  of  the  determination  of  cr  as  well  as  the  interval  of  timesteps 
over  which  values  of  vmax  were  selected. 

Is  there  a  minimum  N  for  which  qeq  requires?  What  is  the  role  of  n?  Can  the  relaxation 
time,  rqeq  be  written  in  terms  of  these  variables?  Are  these  results  based  on  our  attempts  to 
characterize  qeq  using  solid,  liquid,  or  gas  equilibrium  conditions  -  recall  that  granular  materials 
exhibits  phenomenology  of  multiple  states  of  matter.  Does  the  stabilization  of  <  F  >  for  n  >  2.5 
imply  a  sufficiently  good  hard-sphere  model?  It  is  clear  that  much  work  still  needs  to  be  done 
investigating  this  topic. 
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2.4  Mathematical  fit  of  KE  parameter  space 

Coming  back  to  shock  mitigation  capability,  a  mathematical  fit  for  the  KE  =  KE(u),q,N  >  10) 
parameter  space  is  proposed  and  has  the  form: 


EK(u,q)  =  AeBqCeDuE 


(2.16) 


which  corresponds  to  a  two-dimensional  Weibull  distribution.  Due  to  the  gaussian  and  exponential 
natures  of  q  and  lo,  respectively,  the  form  in  equation  (2.16)  then  is  well  suited  when  the  exponents 
are  restricted  to  C  >  1  and  E  <  1.  For  simplicity  we  set  E  =  1  and  C  =  3/2.  The  coefficients 
B  and  D  were  evaluated  using  4th  order  polynomial  fits  and  the  scaling  coefficient  A  is  essentially 
the  point  Ek(N,u>  =  q  =  0).  It  turns  out  that  a  second  order  fit  was  not  sufficiently  robust  and  a 
higher  order  fit  yielded  marginal  gains  at  the  cost  of  mathematical  encumbrance. 

Additionally,  this  model  currently  lacks  the  rigor  sufficient  for  planar  behavior  in  the  limit  of 
small  N.  It  is  unclear  at  this  point  if  N  can  be  completely  decoupled  from  u>  and  q  and  written  as 
an  additional  exponential  modification.  In  all  likelihood,  the  coefficients  B  and  D  would  be  written 
as  functions  of  N. 

In  evaluating  the  fits,  it  was  found  that  KE(lu,  q,N  =  20)  is  described  quite  well  by: 

Ek(uj  q )  =  Q  35544  .  e([-l-5055  10-594+4.01610-493-3.981-10-3g2+0.01479-0.05435]93/2)  _ 

g  ([4. 144-10^594  +  1. 955- 10_3g3-0. 0396292+0. 02887(3-8. 341]w)  ^  47 \ 


It  is  compelling  to  rewrite  (2.17)  in  the  more  suggestive  form: 

(2.18) 
(2.19) 


EK{w,q)  ~ 


eag(n+3)  e0q(n+2)  g79(”+1)  g Sqn  q 


(»— 1) 


where  the  powers  of  q  have  been  adjusted  to  be  written  as  a  series  in  n.  The  occurrence  of  an  apparent 
series  in  n  is  quite  striking.  Note  that  in  (2.17)  for  increasing  strength  of  q,  the  coefficients  decrease 
by  orders  of  magnitude.  These  weaker  weightings  are  mandatory  since  they  are  arguments  in  the 
exponential  and  decay  of  KE  would  occur  too  quickly  were  their  value  too  high.  Additionally,  note 
the  mixed  term  qui  in  (2.17).  This  is  indicative  of  a  many-body  effect:  one  cannot  have  restitution 
in  a  single  particle  chain  since,  by  definition,  it  requires  at  least  two  to  evaluate  the  loading  and 
unloading. 
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Figure  2.26  compares  the  numerical  results  and  the  model  for  the  case  of  N  =  20.  The  results 
are  practically  indistinguishable  except  for  the  minor  timestep  errors  in  the  simulation  which  lead 
to  the  roughness  in  the  difference  plot  for  large  q.  For  this  large  value  of  N  and  q ,  the  tail  particle 
is  rather  small  and  has  a  large  velocity.  The  time  step  is  such  that  kinetic  energy  plots  lose  their 
accuracy  because  the  algorithm  picking  the  1st  peak  is  seeing  less  smoothness  in  the  KE(f)  space. 
As  such,  the  model  provides  a  better  assessment  for  the  state  of  the  system. 


48 


particle 


(a) 


1000 


(b) 


1000 


(C) 


1000 


(d) 


1000 


(e) 


1000 


(f) 


Figure  2.14:  Normalized  kinetic  energy  landscape.  Vertical  axes  represents  the  particle  id  while  the 
horizontal  axis  is  time  out  to  1000  /rs. 
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Figure  2.15:  Velocity  distribution  with  superimposed  Gaussian  fit  for  each  particle  in  a  STC  with 
N  =  4,  q  =  0.  Particle  1  represents  the  end  of  the  chain  and  particle  4  represents  the  head  where 
the  initial  impulse  was  applied. 
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Figure  2.16:  Velocity  distribution  with  superimposed  Gaussian  fit  for  each  particle  in  a  STC  with 
TV  =  20,  q  =  0.  Particle  1  represents  the  end  of  the  chain  and  particle  20  represents  the  head  where 
the  initial  impulse  was  applied. 
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Figure  2.17:  Average  kinetic  energy  of  the  system  as  a  function  of  n  and  N. 
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Figure  2.18:  Average  kinetic  energy  of  the  system  determined  from  simulation  (circles)  and  theoret¬ 
ical  expectation  by  the  Virial  Theorem  (dashed). 
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Figure  2.19:  Partition  of  energy  by  grain  number  for  several  STCs  where  N  =  {3,8,14,20}  and 
q  =  {0,0.05,0.1}  at  t  =  100ns.  The  red  line  indicates  the  theoretical  expectation  that  <  KEn  >= 

n/ (n+2) 
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Figure  2.20:  Partition  of  energy  by  grain  number  for  several  STCs  where  N  =  {3,8,14,20}  and 
q  =  {0,0.05,0.1}  at  t  —  500/its.  The  red  line  indicates  the  theoretical  expectation  that  <  KEn  >= 

n/ (n+2) 
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Figure  2.21:  Partition  of  energy  by  grain  number  for  several  STCs  where  N  =  {3,8,14,20}  and 
q  =  {0,  0.05,  0.1}  at  t  =  1000/xs.  The  red  line  indicates  the  theoretical  expectation  that  <  KEn  >= 
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Figure  2.22:  Fluctuations  as  a  function  of  time  and  n  for  N  —  100,  <7  =  0. 
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Figure  2.23:  (a)  Variation  of  <  F  >  with  N  for  various  n  and  fitting  function,  (b)  Variation  (scatter) 
of  k  with  n.  (c)  Variation  of  A  with  n  and  corresponding  fitting  function,  (d)  Variation  of  B  with 
n  and  corresponding  fitting  function. 
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Figure  2.24:  <  F  >  as  a  function  of  n  for  various  N  determined  from  simulation  (circles).  Derived 
fitting  function  <  F  >=<  F(n,  N)  >  is  superimposed  (dashed)  in  each  panel,  n  =  2  data  is  excluded 
as  it  doesn’t  fit  the  function. 


57 


N=16 


N=18 


N=50 


Figure  2.25:  \vmax\  as  a  function  of  time  for  N  =  16, 18,  20, 50. 
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Figure  2.26:  The  Simulated  and  Modeled  STC. 


Chapter  3 


The  Decorated  Tapered  Chain 
(DTC) 


3.1  Introduction 

The  content  of  this  chapter  was  reported  by  Doney  and  Sen2'’28  in  Physical  Review  Letters.  The 
DTC  (Figure  3.1)  can  be  assembled  from  the  STC  by  introducing  a  single-sized  interstitial  grain  of 
radius,  reven,  between  every  member  of  the  STC  .  We  constrain  the  system  to  an  odd  number  of 
particles  such  that  the  grains  that  form  the  ends  of  the  STC  are  still  the  outer  members  in  the  DTC. 
Additionally,  we  presume  that  these  interstitial  grains  will  be  equal  to  or  smaller  than  the  smallest 
member  of  the  STC  (which  has  radius  rjv);  therefore,  reven  =  /rjv,  where  0  <  /  <  1.0 — although 
the  flexibility  is  already  built  in  for  /  to  be  any  size. 

It  is  immediately  clear  that  the  grain  size  mismatch  changes  as  a  function  of  position  along  the 
DTC.  This  is  in  stark  contrast  to  the  STC  where  successive  grains  are  always  smaller  (or  larger)  by 
the  same  amount.  It  is  possible  then  to  have  DTC  chains  that  appear  to  resemble  monodisperse 
chains  for  only  a  portion  of  the  chain. 

Section  3.2  derives  the  HSA  and  a  fascinating  limiting  case  for  the  DTC.  Normalized  KE  pa¬ 
rameter  spaces  are  also  presented.  These  results  are  followed  in  section  3.3  by  the  numeric  solution 
to  the  EOM  which  includes  a  discussion  on  energy  partitioning.  Supporting  plots  highlighting  the 
instantaneous  energy  breakdown  per  grain  for  a  variety  of  chains  in  10  /is  intervals  is  provided  in 
Appendix  D. 


3.2  Hard-Sphere  Approximation 


In  deriving  an  approximation  for  the  DTC,  the  process  is  more  cumbersome  and  the  conservation 
equations  for  mass  and  energy  are  carried  out  for  several  terms  until  a  pattern  emerges.  Our  primary 
interest  is  in  deriving  an  expression  for  the  normalized  kinetic  energy, 


KEout  _  mN  (V'N\2  _  mN  {  (  v'n  \  ( vi+ 2 \  f  vi+ 1  \ 

KEvn  m1\v1)  TOl  l  \ujv-i/  \vi+i)\  v'i  )  VO// 


(3.1) 
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0.1 
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(a)  f  =  1;  qd  = 


(b)  f  =  1;  qd  = 


(e)  f  -  0.3;  qd  -  0.1  (f )  f  =  0 . 3;  qd  =  0 


Figure  3.1:  Several  example  DTCs  created  by  varying  /  and  qd  for  constant  N  =  13. 

Note  that  tapering  in  the  DTC,  qd,  is  defined  differently  than  for  the  STC.  Thus,  rt+ 2  =  (1  —  qd)ri- 
We  will  eventually  look  for  forms  of  v'l+l/vl  and  then  generalize  for  N  particles  or  (N  —  1)  collisions. 
First,  the  relationship  among  masses  and  radii  must  be  evaluated.  Assembling  the  radii,  we  have 
starting  with  the  largest  bead, 

n 

n+1  =  frN 
Ti+ 2  =  n-  qdTi  =  (1 

n+ 3  =  frN 

^  i-\-  4  =  T  i-\-2  QdT  z+2 

n+ 5  =  frN 

^i+6  =  i+ 4  QdTi-\-  4 

rN-x  =  frN 

rN  =  rN-2  -  qdrN-2  =  (1  -  qd)rN-2  =  e(N~1)/2ri 


-  qd)n  =  ert 


=  (1  -  qd)ri+ 2  =  e2Vi 


=  (1  -  qd)ri+ 4  =  e’ri 
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The  major  equations  for  radii  are  therefore 


rjv  =  e^N  l^2ri 


r(i+l),(i+3),---,(JV- 

-l)  =  /e(W  1)/2n 

Recall  for  masses  that  rm  =  pVj  =  |7r r3p  = 

=  77 rf .  Note  that 

cancel  once  the  conservation  equations  are  put  into  use,  we  will : 

for  nrii  then  become, 

TOi  - 

,  e°r3 

mi+ 1  - 

J  r£+l  —  i 

m-i+2  ~ 

,  e3r3 

m.i+ 3  - 

J  ri+3  —  ^r£ 

mi+ 4  ~ 

,  e6r3 

m,N-i  ~ 

1 

II 

00 

rriN  ~ 

,  f3(JV-l)/2  3 

i 

(3.2) 


(3.3) 


where  A  =  y3e3(JV-i.)/2.  These  relations  can  now  be  used  to  set  up  the  conservation  equations. 
Beginning  with  momentum  and  assuming  that  each  subsequent  particle  in  the  chain  begins  at  rest, 
one  can  solve  for  the  first  five  collisions.  Primes  and  double-primes  indicate  post-collision  states. 
A  primed  quantity  denotes  the  first  post-collision  state  of  a  sphere.  That  subsequent  sphere  will 
serve  as  the  input  to  the  next  collision.  To  keep  track  of  its  velocity  after  the  second  collision,  it  is 
denoted  by  a  double-prime  and  will  eventually  be  eliminated.  With  e  =  (1  —  qd), 


rmvi  =  mWi  +  mi+1v'i+1 
l^i+1  =  raz+1^+1  “1“  ^£+2^+2 
^i+2^+2  =  mi+‘2vi+2  +  mi+Svi-\-3 
^ i+3vi-\-3  =  ^i+3^+3 
^z+4^-|-4  =  H”  ^£+5^+5 


Vi  =  Vj  +  Avm 
M+i  =  Av’{+1  +  e3u'+2 
eau(+2  =  ea<+2  +  Av^+a 
AVj+ 3  =  Av”+  3  +  ebu'+4 
ebu(+4  -  ebv”+4  +  Av’i+b 


(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 
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From  the  pattern  in  equations  (3. 4-3. 8),  (3.4)  can  be  rewritten  as  e°v'i  =  e°v'i  +  Av'i+1.  An  evaluation 
of  energy  conservation  yields  the  same  form  as  equations  (3. 4-3. 8)  except  velocities  are  squared: 


v2=v'2  +  Av'2+l 

(3.9) 

Av?.  (  =  ^4+i  +  e3v?+2 

(3.10) 

,3/2  —  f3v"2  -1-  At/2 

e  vi+ 2  —  e  vi+ 2  +  Avi+ 3 

(3.11) 

4-/2  _  4- "2  ,  6  /2 

Avi+ 3  —  Avi+3  +  e  vi+ 4 

(3.12) 

,6,/2  —  r61>"2  -1-  4A2 

e  A+4  —  e  vi+ 4  +  ^A+5 

(3.13) 

Next,  combine  equations  (3. 4-3. 8)  and  (3.9-3.13)  to  eliminate  the  double-primed  terms  and  form 
the  velocity  ratios:  ^4 — ,  etc.  Beginning  with  (3.4),  isolate  v\  and  square  it  to  obtain 

v’i  =  vf  —  2 Aviv'i+1  +  a2v'2+1.  Next,  substitute  this  into  (3.9)  and  rearrange  to  obtain  This  is 
then  repeated  for  equations  (3.5,3.10),  (3.6,3.11),  etc.  to  obtain  the  following  ratios: 


v 


*+ 1 


Vi 


V 


V 


%- (-2 
7 

i+ 1 
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i+3 


Vi+ 2 


V 


V 


i-\-  4 

7 

i+ 3 


2e° 

(3.14) 

A  +  e° 

2A 

(3.15) 

e3  +  A 

2e3 

(3.16) 

A  +  e3 

24 

e6  +  A 

(3.17) 

where  hindsight  has  allowed  us  to  insert  terms  of  e°  in  (3.14).  With  an  eye  fixed  on  (3.1),  results 
may  be  merged.  Thus 
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2  A 


2e3  \  (  2 A  \  (  2e6 
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2  A 


(3.19) 


N= 3 


N=5 


N=7 


The  ratio  can  be  put  into  closed  form  to  obtain, 


UN 

Vi 


{N^(  2e>W-1) 


2A 


l=i 


A  +  e3(J_1)  /  \  e3l  +  A 


(3.20) 


(JV-l)/2 

2JV-1  4(JV-l)/2  TT  |  _ _ 

U  +  £3(i-i) 
1=1  \ 


1 


e3i  +  A 


(3.21) 


Turning  to  the  mass  ratios  and  using  expressions  from  (3.3),  it  appears  that  most  terms  cancel: 


mN 

mi 


leading  to  the  simple  expression, 


(3.22) 

(3.23) 


^  =  e3(iV-i)/2  (3.24) 

mi 


One  can  now  identify  the  normalized  kinetic  energy  (3.1)  by  squaring  (3.21)  and  combining  it 
with  (3.24)  to  form, 
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Figure  3.2:  Normalized  kinetic  energy  surfaces,  Ek  =  Kout/Kin,  for  the  decorated  chain  under 
the  hard  sphere  approximation  as  functions  of  the  number  of  spheres,  N,  fractional  size  of  interstitial 
sphere,  /,  and  tapering,  qd- 


rp(dtc) 


(N- 1)/2 


(4.4f3/2)^v  n 


P3(t-1) 


L=A  (A  +  e30'-1))(e3i  +  A) 


(3.25) 


These  results  are  plotted  in  Figure  3.2.  The  effects  of  the  interstitial  sphere  are  remarkable  when 
compared  to  the  simple  chains  in  figure  2.2(a).  For  a  modest  value,  /  =  0.7,  it  takes  very  few 
particles  to  reduce  the  outgoing  kinetic  energy  considerably.  E^tc>  decays  as  a  gaussian  or  sigmoid 
with  increasing  qd,  and  exponentially  with  increasing  N. 

It  is  difficult  to  draw  any  physical  intuition  from  (3.25).  However,  a  very  curious  and  astonishing 
result  occurs  in  the  limit  qd  =  0: 

ekU)\„-0  =  (|75^p)  ■  (3-26) 


This  limit  is  equivalent  to  equation  (2.7)  under  the  exchange  /  <£=>  (1  —  qs). 

It  is  clear  that  /  =  1  should  imply  qs  =  0  since  they  both  generate  monodisperse  chains.  That 
this  equivalency  goes  beyond  that  special  case  is  quite  unexpected.  One  can  now  begin  to  see  the 
incredible  effect  /  has  on  the  energy  mitigation  capability  when  an  infinite  potential  is  invoked:  for 
/  =  0.3  —  a  typical  value  we  might  consider  —  the  equivalent  tapering  in  the  simple  chain  would  be 
qs  =  0.7.  This  value  is  7  times  larger  than  any  system  we  had  previously  considered  and  could  be  a 
significant  system  integration  challenge.  Visually,  for  hard-spheres,  the  energy  mitigation  capability 
of  the  simple  chain  shown  in  Figure  2.1  ( qs  =  0.1)  is  identical  to  that  for  a  decorated  chain  similar 
to  that  shown  in  Figure  3.1(b)  but  with  qd  =  0,  N  =  10,  /  =  0.9. 
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3.3  Numeric  Solution 


3.3.1  KE  Parameter  Space  Behavior 

Figure  3.3  highlights  the  computational  results  for  the  decorated  chain.  Recall  that  the  inertial 
mismatch  between  neighboring  grains  in  decorated  chains  change  as  a  function  of  position  along 
the  chain.  This  is  what  is  believed  to  be  the  cause  of  a  ripple  in  the  surface  of  the  Ek  plots  that 
propagate  toward  the  origin  as  /  decreases.  As  one  might  expect,  such  behavior  is  a  function  of 
N,  qd ,  and  /.  The  effect  vanishes  for  /  <  0.6,  approximately.  At  about  this  threshold,  the  interstitial 
grain  is  not  much  smaller  (less  massive)  than  the  grains  toward  the  end  of  the  chain.  The  explanation 
is  that  as  an  impulse  propagates,  energy  transmission  becomes  increasingly  efficient  due  to  smaller 
inertial  mismatches  —  a  prerequisite  for  admitting  solitary  waves.  Thus  the  system  changes  from  a 
shock  absorber  to  shock  transmitter.  This  effect  however  must  compete  with  compressive  effects  in 
some  manner  since  no  such  behavior  is  present  for  hard  spheres  even  though  it  too  has  a  position- 
dependent  inertial  mismatch. 

The  hard  sphere  approximation  grossly  exaggerates  the  shock  mitigation  capability  of  the  deco¬ 
rated  chain.  Additionally,  it  doesn’t  pick  up  the  surface  feature  resulting  from  a  competition  between 
particle  overlap  nonlinearity  and  variable  inertial  mismatches  between  neighboring  grains.  Thus  a 
hard-sphere  analysis  is  inadequate  for  the  DTC.  Simulations  suggest  that  for  /  =  0.3,  N  =  5,qd  =  0.1 
(panel  f),  one  can  disperse  energy  within  the  chain  such  that  only  about  10%  of  that  put  into  the 
system  is  transmitted  to  the  end  with  the  initial  pulse. 

3.3.2  Energy  Partition 

A  simple  view  of  energy  partitioning  in  the  DTC  is  not  possible  given  the  vast  number  of  possible 
chain  configurations,  belittling  those  of  the  STC.  Changing  N  has  a  much  more  severe  impact  on 
the  results  because,  by  design,  it  affects  the  results  everywhere  in  the  chain.  For  example,  if  we  take 
the  mass  ratio  of  interstitial  grains  for  N  =  21  and  N  =  11  chains,  with  qd  and  /  identical,  the 
results  scale  astonishingly:  mjv=2i/wjv=n  oc  (1  —  qd)15'  For  large  tapering,  qd  =  0.1,  interstitial 
grains  in  a  chain  with  N  =  21  have  about  one- fifth  the  mass  as  their  counterparts  where  N  =  11. 

Rather  than  focusing  on  the  system,  it  is  of  interest  to  probe  how  a  DTC  partitions  the  energy 
among  particles  as  an  impulse  propagates.  Figure  3.4  illustrates  the  instantaneous  kinetic  energy 
per  grain  for  various  configurations  at  t  =  52fis  where  qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
Also  included  in  each  panel  is  a  silhouette  of  the  specific  chain  and  the  total  kinetic  energy  of  the 
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Figure  3.3:  Numerically-produced  normalized  kinetic  energy  surfaces,  Ek  =  Kout / Kin ,  for  the 
decorated  chain  as  functions  of  the  number  of  spheres,  N,  fractional  size  of  interstitial  sphere,  /, 
and  tapering,  qj.  Several  sample  chains  are  identified  in  panels  (d-i). 

system  as  a  function  of  time.  The  black  disk  indicates  the  current  simulation  time.  It  should  be 
noted  that  panel  (a)  is  identical  to  figure  2.13(a)  since  they  are  both  monodisperse  chains.  A  cursory 
glance  of  all  panels  reveal  that  the  effect  of  qd  is  to  spread  the  impulse  out  over  many  grains  in  a 
manner  similar  to  the  STC.  And  secondly,  when  decreasing  /  the  energy  appears  to  be  distributed 
to  the  larger,  non-interstitial  (odd  numbered)  grains.  It  is  hypothesized  that  additional,  neighboring 
interstitial  spheres  would  further  separate  the  energy  along  the  chain.  What  this  also  appears  to  do 
is  turn  the  DTC  into,  effectively,  a  binary  collision  system  since  the  amplitudes  of  interstitial  spheres 
is  quite  superficial.  Also  of  note  is  that  the  speed  of  energy  transmission  appears  more  dependent 
on  qd  than  /. 

The  division  of  energy  into  the  larger  grains  for  small  /  appears  to  be  a  result  of  their  larger 
masses,  rather  than  the  increased  “rattling”  of  interstitial  spheres.  For  example,  the  mass  ratio 
between  grains  15  and  14  in  panel  (i)  is  about  0.004.  This  value  becomes  more  matched  as  one 
approaches  grains  1  and  2  which  has  a  ratio  of  about  0.03.  Even  though  the  even-numbered  grains 
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Figure  3.4:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  52/us  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 


have  a  much  higher  velocity  than  their  larger  neighbors,  kinetic  energy  only  scales  as  v2  versus 
m  ~  r3 .  The  smaller  system  kinetic  energy  plots  in  each  panel  reveals  the  complicated  nature  of  the 
system — a  consequence  of  the  competition  among  f,qd,N.  The  dynamics  tend  to  be  more  smooth 
for  small  /  because  of  the  smaller  role  played  by  interstitial  spheres  to  the  total  system  energy.  One 
area  that  still  needs  investigating  is  possibility  of  multiple  grains  overlap  for  very  small  /.  In  this 
case,  the  EOM  would  need  to  be  modified  to  take  into  account  the  next  layer  of  neighbors. 

Panels  (b,c)  may  not  necessarily  admit  SWs,  but  localized  energy  propagation  does  occur.  In 
both  cases,  the  amplitude  dampens,  in  agreement  with  Manciu59 — this  is  more  obvious  with  larger 
N  (not  pictured).  Panels  (f,i)  quickly  spread  out  the  energy  because  of  the  finite  tapering  which  is 
also  the  reason  for  an  increasing  KE  ramp  (shown  in  the  inset).  However,  tapering  qd  now  must 
compete  with  /  as  results  vary  significantly  among  other  panels.  A  thorough  analysis  of  that  rivalry 
has  not  yet  been  investigated. 

It  is  instructive  to  observe  in  time,  how  waves  propagate  through  such  systems.  Appendix  D 
shows  a  sequence  of  these  plots  in  10  /rs  intervals.  One  can  flip  through  the  pages  so  that  they 
appear  animated. 
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Chapter  4 

Hydrocode  Simulations  of  Tapered 
Chains 


There  are  two  reasons  that  stand  out  when  considering  the  simulation  of  TCs  in  a  hydrodynamics 
code  (hydrocode).  The  first  is  in  probing  the  regime  where  plastic,  and  therefore  permanent,  de¬ 
formation  can  occur.  Secondly,  it  represents  the  next  step  towards  the  development  of  a  potential 
TC  armor  panel.  This  brief  chapter  describes  preliminary  efforts  in  running  such  2D  and  3D  simu¬ 
lations.  These  results  were  presented  at  the  March  2006  meeting  of  the  American  Physical  Society 
in  Baltimore,  MD. 

ALEGRA 18,79,110 — or  Arbitrary  Lagrangian  and  Eulerian  General  Research  Application — is  a 
massively  parallel,  finite  element  hydrodynamics  code  with  modular  physics  capability  developed 
by  Sandia  National  Laboratories.  The  general  procedure  for  running  such  simulations  consists  of 
mesh  generation,  simulation,  and  post-processing.  The  first  step  is  performed  using  CUBIT22  which 
discretizes  the  computational  space.  At  each  node,  the  conservation  equations  are  solved  using 
ALEGRA.  Finally,  visualization  is  performed  using  Ensight17.  While  the  manuals  of  the  code  are 
publicly  available,  the  software  itself  is  export-controlled. 

Both  lagrangian  and  eulerian  meshes  have  been  considered  for  3D  and  2D  meshes,  respectively. 
In  the  latter  case,  material  advects  through  a  stationary  mesh1.  This  is  in  contrast  to  the  Lagrangian 
framework  where  the  material  is  the  mesh.  An  inherent  difficulty  in  running  Eulerian  calculations, 
is  that  mixed  cells  are  endemic  to  the  mesh.  That  is,  any  cell  can  be  occupied  by  more  than  1 
material:  part  air,  part  water  for  example.  As  time  progresses  and  the  material  advects  through 
the  mesh,  a  code  needs  to  keep  track  of  that  boundary  between  materials  which  has  to  be  deter¬ 
mined  at  a  resolution  finer  than  the  mesh.  This  is  known  as  interface  tracking  (IT)  or  interface 
reconstruction  (IR).  Currently,  ALEGRA  supports16  three  schemes:  SLIC2,  SMYRA38,  and  NEW 
SMYRA.  Results  from  simulations  should  be  independent  of  the  mesh  ( “mesh  convergence” )  as  well 
as  the  IT  scheme,  much  in  the  way  that  physics  is  independent  of  a  coordinate  system.  In  Lagrange 

1  More  specifically,  there  is  a  lagrange  step  where  the  mesh  deforms  to  follow  the  material  and  then  is  mapped  back 
into  the  eulerian  mesh 

2Simple  Line  Interface  Calculation 
3  Sandia  Modified  Young’s  Reconstruction  Algorithm 
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Figure  4.1:  Computational  domain  and  problem  setup  for  the  ALEGRA  simulation. 


frameworks,  there  are  no  mixed  cells  —  so  it  is  the  preferred  method  assuming  deformations  are  not 
too  large  to  cause  element  inversions. 

Properties  of  inert  and  reactive  materials  and  equations  of  state  have  been  evaluated  semi- 
empirically  for  the  past  half-century.  Many  of  those  results  are  available  as  tabular  databases  and 
can  be  read  by  various  computer  codes.  The  SESAME  EOS  is  one  such  behemoth  and  used  frequently 
in  hydrocode  simulations.  One  also  needs  a  material  strength  model  and  the  Johnson-Cook16  model 
is  common, 


Y=  [A  +  BeN] [1  +  C In  (max(0.002,  e))] 


T  —  Tr 
T  —  T 

m  r 


m 


where  the  A,  B,  C,  N,  to,  T  parameters  are  material  dependent,  e  is  the  equivalent  plastic  strain, 
and  Tr,Tm  are  the  room  and  melt  temperatures,  respectively.  The  values  of  these  constants  have 
been  evaluated  for  a  great  many  materials.  Note  that  all  material  models  are  not  created  equal. 
Some  models  perform  well  through  phase  transitions  for  example,  while  others  may  not,  i.e.  Mie- 
Griineisen. 


4.1  2D  Eulerian  ALEGRA  simulations 

The  initial  effort  consisted  of  inserting  copper  disks  into  a  2D  eulerian  mesh  with  cartesian  symme¬ 
try.  This  implies  that  the  problem  is  really  that  of  infinitely  long  cylinders  since  the  symmetry  is 
translational  along  the  .z-axis.  Copper  disks  were  used  because  the  elastic-plastic  behavior  in  the 
SESAME  model  is  well-characterized  (even  across  multiple  phases)  and  has  been  utilized  for  quite 
some  time.  Figure  4.1  highlights  the  problem  setup.  The  geometry  consists  of  a  5-particle  monodis- 
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Figure  4.2:  Evolution  of  \crxx\.  Color  scale  is  logarithmic. 


perse  chain  adjacent  to  a  large  wall.  Grain  1  is  impacted  by  a  smaller  sphere  moving  to  the  right 
with  velocity,  v.  The  purpose  of  LTi  is  to  measure  the  grain-grain  interaction  and  obtain  a  value  for 
Sn  while  LT2  can  record  force.  The  incorporation  of  such  Lagrange  tracer  particles  —  which  move 
with  the  material  —  allows  the  user  to  track  various  quantities  over  time  without  worrying  about 
displacement.  This  problem  is  ideal  for  a  first  series  of  verification  and  validation  studies  that  could 
be  based  on  the  work  of  Rossmanith  and  Shukla91.  The  following  steps  would  be  required:  increase 
the  number  of  spheres,  add  a  gravitational  field  and  change  the  material  parameters.  For  the  case 
of  increasing  zig-zag  obliquity,  the  disks  would  need  to  be  put  between  walls  whose  distance  changes 
with  the  angle  of  obliquity. 

Preliminary  results  of  a  test  simulation,  where  \crxx\  —  the  absolute  normal  component  of  stress 
along  x  in  units  of  Pascals,  are  shown  in  Figure  4.2.  The  input  velocity  is  about  2.5  mm//zs  and 
the  color  scale  is  logarithmic.  The  regions  in  the  top  3  panels  —  in  what  looks  like  ringing  —  are 
believed  to  be  artifacts  of  an  improper  lower  boundary  in  the  color  map.  It  is  therefore  most  likely 
numerical  noise.  What  is  clear  however,  is  that  the  shock  propagation  adds  a  level  of  fidelity  not 
available  in  the  analysis  of  earlier  chapters. 
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Figure  4.3:  Basic  3D  simulation  of  a  TC.  Note  that  the  bulk  behavior  of  the  spheres  is  consistent 
with  expectation  even  though  dynamics  take  place  on  a  elemental  or  nodal  basis.  A  small  amount 
of  plastic  flow  is  visible. 

4.2  3D  Lagrangian  ALEGRA  simulations 

More  interesting  to  the  research  and  where  the  effort  is  now  focused,  is  the  simulation  of  a  3D 
alignment  of  spheres.  The  jump  in  complexity  however  is  quite  considerable  and  requires  the  incor¬ 
poration  of  a  global  contact  model.  The  contact  mechanics  package  ALEGRA  utilizes  is  known  as 
ACME13.  This  area  is  probably  the  most  critical  in  setting  up  the  problem  and  therefore  the  one 
that  will  require  the  most  care. 

In  the  3D  environment,  individual  spherical  meshes  are  created  and  positioned  similar  to  how 
a  tapered  chain  is  parameterized.  An  initial  simulation  was  run  to  verify  that  sphere  velocities 
increase  along  a  chain  when  finite  tapering  is  introduced.  That  result  is  presented  in  figure  4.3 
where  N  =  4,  q  ~  0.05.  In  the  animation,  it  is  clear  that  the  velocity  of  the  right-most  particle  is 
larger  than  the  trailing  spheres.  The  same  is  true  for  the  penultimate  sphere,  etc.  Also  visible  is  a 
decreasing  amount  of  deformation  in  the  spheres  as  one  moves  rightward.  In  this  sense,  the  spheres 
have  absorbed  energy  even  though  plastic  flow  is  evident.  What’s  important  to  point  out  is  that 
even  though  interactions  are  being  performed  on  a  nodal  or  elemental  basis  where  the  latter  is  of 
arbitrary  size,  the  bulk  behavior  of  the  system  is  consistent  with  expectations.  Upon  close  inspection 
one  can  gauge  the  resolution  of  the  mesh.  For  these  problems,  there  are  50  azimuthal  divisions  and  6 
radial  divisions  —  corresponding  to  1200  elements  —  per  particle.  The  CUBIT  script  which  creates 
the  spherical  meshes  is  listed  in  appendix  F  while  the  ALEGRA  input  script,  which  specifies  the 
relevant  physics  is  listed  in  appendix  G.  Note  that  in  both  cases  the  real  input  for  both  CUBIT 
and  ALEGRA  is  much  larger.  Both  scripts  have  utilized  a  package  known  as  APREPRO  which  is  a 
preprocessing  utility  that  allows  one  to  incorporate  algebraic  expressions. 

It  was  useful  to  quantitatively  compare  the  results  between  numerically  solving  the  EOM  from 
chapter  2  and  ALEGRA  solutions.  An  exact  match  was  not  expected  since  the  latter  has  barely 
been  tested  and  refined.  Figure  4.4  reflects  this.  Both  simulations  consisted  of  a  TC  where 
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N=5 ,  q=0 . 07 ,  frictionless 


Figure  4.4:  Velocity  of  the  head  particle  as  a  function  of  time  for  both  EOM  and  ALEGRA  simula¬ 
tions  where  N  =  5,  q  =  0.07  and  restitutive  losses  are  neglected. 

N  =  5 ,q  =  0.07.  In  the  case  of  the  EOM  problem,  restitutive  losses  were  neglected.  None  were 
specifically  inserted  into  the  contact  model  for  the  3D  simulation;  however,  it  is  unclear  whether 
there  is  a  nonzero  default  value.  In  addition,  the  boundary  conditions  are  not  the  same  as  the  EOM 
problem.  In  fact,  no  boundary  conditions  along  the  direction  of  motion  were  specified,  so  the  parti¬ 
cles  can  continue  indefinitely  —  which  is  why  the  ALEGRA  curve  shows  a  saturation  that  continues 
unabated.  Also,  it  is  unclear  why  there  is  a  step-like  decay.  It  is  very  possibly  that  mesh-resolution 
is  playing  a  role  where  the  number  of  radial  elements  would  have  the  biggest  impact.  Another 
point  is  that  the  material  models  are  different  between  the  two  cases.  While  many  metals  can  be 
approximated  by  a  Poisson  ratio  of  0.33,  their  strengths  can  vary  substantially. 

As  an  example  of  the  capabilities  of  the  code  and  the  last  simulation  performed  to  date,  the 
initial  velocity  was  substantially  increased  to  500  m/s.  Snapshots  of  the  results  are  visible  in  figure 
4.5.  In  these  panels,  color  is  scaled  linearly  by  velocity  with  dark  blue  representing  spheres  at 
rest.  The  amount  of  deformation  is  significant  and  permanent.  Again,  these  spheres  are  copper  and 
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simulations  have  not  yet  been  performed  with  other  metals.  In  the  last  panel,  the  copper  “pancakes” 
are  still  moving  at  approximately  250  m/s. 

As  a  final  note,  many  interesting  questions  can  now  be  asked  and  explored  with  the  use  of 
ALEGRA.  While  testing  of  the  systems  and  proper  and  careful  setup  still  need  to  be  performed,  there 
is  now  the  capability  to  investigate  the  combined  effects  of  many  TCs.  Antiparallel  arrangements 
can  be  inserted  into  the  mesh  and  surrounded  by  other  metals,  insulators,  etc.  These  can  be  encased 
by  various  alloys  of  steel  and  struck  with  a  flyer  plate  (see  figure  5.1.  The  resulting  force  can  be 
measured  on  the  output  by  surface  tools  in  the  visualization  software.  In  addition,  (virtual)  tracer 
particles  can  be  put  at  arbitrary  locations  to  monitor  local  effects.  The  possibilities  are  actually 
quite  staggering  —  stay  tuned. 
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Figure  4.5:  Time  elapsed  sequence  where  the  initial  velocity  of  the  left-most  grain  has  been  increased 
to  500  m/s  (red).  The  remaining  spheres  are  initially  stationary  (blue). 
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Chapter  5 

Closing  Remarks 


On  March  12,  2006,  Joe  Ray  —  a  close  family  friend,  employee  of  Buncombe  County  in  Asheville, 
NC,  and  an  explosives  ordnance  specialist  —  was  clearing  mines  for  civilian  and  coalition  traffic 
in  Asadabad,  Afghanistan  during  Operation  Enduring  Freedom.  The  up-armored  humvee  that 
was  carrying  him,  Sgt.  Kevin  Akins,  Spc.  Joshua  Hill,  and  Sgt.  Anton  Hiett  of  the  391st  Engineer 
Battalion4,12  was  hit  by  a  stationary  and  camoflauged  improvised  explosive  device.  No  one  survived. 
Such  vehicles  are  overmatched  by  the  seemingly  endless  quantity  of  explosive  that  can  be  put  into 
such  contemptible  devices.  It  is  just  one  example  —  and  a  personal  one  —  where  there  will  always 
be  a  need  to  improve  our  ability  in  absorbing  blast  energy.  While  the  technology  that  has  been 
discussed  in  this  dissertation  is  not  claiming  to  be  able  to  mitigate  the  meganewtons  of  force  that 
such  blasts  can  generate,  it  represents  the  basic  research  which  hopefully  paves  the  way  for  us  to  get 
there.  To  that  end,  we  have  investigated  the  shock  mitigation  and  nonlinear  dynamical  properties 
-  as  it  pertains  to  the  partition  of  energy  —  of  tapered  metal  spheres  arranged  in  one-dimensional 
arrays. 

The  analysis  has  been  carried  out  where  two  chain  geometries  are  of  particular  interest:  the 
simple  tapered  chain  (STC:  figure  2.1)  and  the  decorated  tapered  chain  (DTC:  figure  3.1).  These 
systems  exhibit  strongly  nonlinear  behavior  due  to  the  Hertz  potential  which  has  been  used  to  govern 
the  contact  mechanics  in  solving  the  equations  of  motion  for  each  sphere.  The  dynamics  and  shock 
mitigation  properties  of  each  chain  have  been  analyzed  as  functions  of  the  relevant  parameters. 
These  are  constant  for  each  chain  and  are  defined  as  the  number  of  spheres,  the  tapering,  and  the 
coefficient  of  restitution.  For  the  DTC,  there  is  an  additional  parameter,  /,  which  governs  the  size  of 
interstitial  spheres  used  to  increase  the  inertial  mismatch  between  neighboring  grains.  The  following 
conclusions  with  a  short  summary  is  offered. 
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Both  chain  geometries  decimate  impulses  substantially. 


It  is  clear  from  the  experiments  performed  by  others  and  the  simulations  and  data  presented  here 
that  both  the  STC  and  DTC  decimate  impulses.  Further,  the  DTC  can  outperform  the  STC  in  fewer 
spheres,  making  it  desirable  in  applications  where  little  space  is  available.  The  primary  mechanism 
responsible  is  symmetry  breaking,  or  inertial  mismatches  between  neighboring  grains.  Essentially, 
momentum  and  energy  remain  trapped  in  the  donor  sphere  thus  spreading  an  impulse  out  in  time  and 
space  as  it  propagates  through  a  chain.  This  is  repeated  at  each  interface  where  there  is  a  mismatch. 
For  equally-sized  spheres,  solitary  waves  (SW)  propagate  and  symmetry  is  maintained  except  at  the 
boundary.  Once  the  boundary  is  reached,  the  SW  is  reconstructed,  but  imperfectly.  That  left-over 
energy  goes  into  the  creation  of  secondary  solitary  waves  (SSW)  which  have  the  same  width,  but 
travel  at  smaller  velocities  due  to  the  small  amplitudes.  This  process  continues  indefinitely.  In  the 
DTC,  the  effect  of  qd  is  to  spread  the  impulse  out  over  many  grains  in  a  manner  similar  to  the  STC. 
Secondly,  the  role  of  decreasing  /  is  it  distribute  energy  to  the  larger,  non-interstitial  grains.  By 
increasing  the  number  of  interstitial  grains,  it  is  hypothesized  that  energy  would  continue  to  spread 
out  and  remain  trapped  in  the  larger  grains.  This  has  applications  for  redirecting  the  energy  if  those 
sites  were  sufficiently  coupled  to  some  transit  channel  —  making  the  problem  2D  or  3D.  One  could 
also  potentially  increase  dissipation  with  clever  selection  of  materials. 

A  hard-sphere  approximation  is  valid  for  the  STC  but  invalid  for  the  DTC. 

Normalized  energy  and  force  surfaces  were  evaluated  to  measure  each  chain’s  trapping  ability.  In 
considering  the  STC,  results  indicate  that  there  is  a  sigmoidal  or  gaussian  decay  on  q  and  an 
exponential  dependence  on  N.  This  trend  is  consistent  with  the  results  from  numerically  solving 
the  equations  of  motion.  For  the  DTC  however,  the  hard-sphere  analysis  and  numerical  solution  are 
markedly  different  in  terms  of  the  former  missing  phenomenology  and  greatly  over-predicting  energy 
absorption. 

The  DTC  hard-sphere  approximation  reduces  to  the  STC  in  a  special  limit. 

It  was  quite  surprising  that  in  the  limit  qd  =  0,  equation  (3.25)  —  which  describes  the  shock 
absorption  quality  of  a  DTC  for  N,  q,  f  —  reduces  to  equation  (2.7)  under  the  exchange  /  *£=> 
(1  —  qs).  This  says  that  a  hard-sphere  chain,  consisting  of  an  alternating  series  of  radii  (where 
'r small  =  ffiarge ),  has  the  kinetic  energy  absorption  equivalency  of  a  STC  with  tapering  qs. 
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Figure  5.1:  Schematic  of  what  a  tapered  chain  armor  panel  might  look  like. 

A  tapered  chain  “armor  panel”  can  be  designed  which  consists  of  anti-parallel  arrays 
of  chains. 

Results  indicate  that  impulses  are  attenuated  regardless  of  the  direction  of  propagation  in  a  STC, 
therefore  one  could  design  a  tapered  chain  “armor  panel”  by  placing  anti-parallel  arrays  within  a 
supporting  matrix.  In  that  sense,  the  collective  abilities  of  hundreds  of  such  TCs  could  be  combined. 
A  schematic  is  outlined  in  figure  5.1.  The  quantities  illustrated  were  used  in  order  to  obtain  an 
estimate 26  of  the  specific  absorbed  energy  of  the  panel  and  yielded  results  that  are  on  the  order  of 
other  technologies  (i.e.  up  to  80  J/g).  The  value  for  the  mass  of  the  panel,  M,  however  is  highly 
subjective  and  the  research  needs  to  proceed  to  a  more  advanced  experimental  state  in  order  to 
refine  the  calculation. 

Initial  3D  simulations  of  STCs  hint  at  continued  energy  absorption  well  beyond 
purely  elastic  behavior. 

Initial  3D,  lagrangian  simulations  have  begun  using  the  modern  and  “massively”  parallel  hydrody¬ 
namics  code,  ALEGRA,  developed  by  Sandia  National  Laboratories.  Results  suggest  that  while  an 
increased  impulse  leads  to  plastic  deformation,  there  are  indications  of  continued  energy  absorption. 
Verification  and  validation  studies  as  well  as  calibration  of  the  global  contact  model  needs  to  be 
performed. 
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An  expression  for  the  period  of  a  single  particle  in  the  overlap  potential  has  been 
evaluated. 

Equation  (1.11)  quantifies  the  period  of  a  single  particle  in  the  overlap  potential,  of  arbitrary  power 
n  and  input  velocity,  ry. 

STCs  do  not  equipartition  energy  but  they  do  satisfy  the  virial  theorem. 

Equipartition  can  be  spoken  of  either  at  a  system  level  or  for  individual  particles.  In  the  former, 
equipartition  is  not  satisfied,  but  discrete  systems  with  arbitrary  n  in  the  overlap  potential  were 
found  to  obey  the  virial  theorem.  At  a  particle  level,  certain  chains  showed  a  tendency  to  equally 
partition  the  available  KE  to  all  members  after  sufficient  time  had  transpired. 

Even  after  long  simulation  time,  STCs  still  fluctuate  about  the  system’s  mean  kinetic 
energy.  They  appear  to  settle  into  a  quasi-equilibrium. 

While  the  mean  value  of  KE  in  STCs  was  found  to  obey  the  virial  theorem  and  particle  veloci¬ 
ties  demonstrate  a  Gaussian  profile,  the  system  settles  into  a  volatile  state  referred  to  as  quasi¬ 
equilibrium.  The  fluctuations  remain  quite  measurable  even  over  long  simulation  time  and  their 
mean  value  has  been  quantified  as  a  function  of  the  number  of  particles  and  the  exponent  in  the 
overlap  potential  (equation  2.13). 

A  great  number  of  questions  still  need  to  be  investigated. 

In  concluding  this  dissertation,  I  pose  the  following  research  areas  for  further  study  beyond  the 
endless  topics  and  questions  that  the  3D  simulations  pose. 

•  A  close  inspection  of  figure  2.23(a)  for  n  =  2,  2.5  indicates  that  equation  (2.13)  governing  the 
decay  in  the  mean  fluctuations  may  not  be  sufficiently  robust.  It  would  be  useful  to  run  the 
simulations  where  N  =  500  or  N  =  103  to  improve  the  formulation.  Also,  why  does  <  F  > 
appear  to  be  independent  of  ql 

•  Daraio23  has  experimentally  verified  shock  absorption  in  chains  where  the  material  properties 
have  been  varied.  What  can  be  learned  from  numerical  simulations,  say  for  material  and  size 
variation? 

•  Is  it  possible  to  exploit  the  phenomenon  of  inelastic  collapse  for  shock  mitigation  purposes? 
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•  What  is  the  proper  way  to  characterize  quasi-equilibrium?  Should  it  be  based  on  equilibrium 
conditions  of  a  particular  state  of  matter  or  combinations  thereof? 

•  What  does  the  stress-strain  curve  of  a  TC  look  like? 

•  What  is  the  best  filler  material  in  a  tapered  chain  armor  panel? 
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Appendix  A 


STC  Code  from  Pfannes80 


/*  PROGRAM  taperch.ain.cpp 


ORIGINAL  VERSION:  J.M.M.  Pfannes 

This  program  consideres  an  one  dimensional  chain  of  spheres  that 
shrink  succesively  in  radius  ("tapered  chain").  Initially  the  spheres 
are  barely  in  contact,  i.  e.  they  just  touch  each  other  and  are  not 
compressed  (zero  loading) . 

The  chain  ends  at  both  edges  at  fixed  walls. 

The  program  calculates  the  interaction  of  the  system  once  disturbed 
by  an  instantaneous  (delta)  impulse  exerted  on  one  end  of  the  chain. 
Restitution  both  between  the  spheres  and  between  the  edge  spheres 
and  the  corresponding  wall  can  be  introduced. 

The  program  does  not  consider  gravity. 

The  EOM  are  solved  with  the  Velocity  Verlet  algorithm. 

scale  of  problem:  mm-mg-musec  (mimimu) 
in  this  scale  unit  of  force:  1000  N 
in  this  scale  unit  of  energy:  1  J 

*/ 


#include  <cmath> 
#include  <iostream> 
#include  <fstream> 
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#include  <cstdlib> 


#include  <string> 

#include  <sstream> 

using  namespace  std; 

/**************************  ALTERABLE  PARAMETER:  ****************************/ 
const  int  nptles=20;  //  total  number  of  particles 

const  double  rh.o=3 . 2  /*  SiC  (mg/mm~3)  */ ,  D=0 . 00326603139013  /*  (mnT2/N)  */ ; 

const  double  rlarge  =5.0;  //  (radius  of  large  ptle  (mm)) 

const  double  q  =  0.0;  //  (tapering  factor  (%)) 

const  double  xn  =  2.5;  //  (exponent  in  potential) 

const  double  dt  =  0.00001;  //  (timestepwidth.  (musec)) 

const  unsigned  int  nsteps  =  100000000;  //  (#  steps  integration  loop) 

const  int  idiagp  =  20000;  //  (stepwidtb  diagnostics) 

const  int  idump  =  20000;  //  (stepwidth.  dump) 

const  double  vlin  =0.0;  //  (initial  v  small  ptle  (mm/musec)) 

const  double  vnin  =  -0.01;  //  initial  v  large  ptle  (mm/musec)) 

const  double  epsilon  =1.0;  //  ((1-  restitution  factor)  all  ptles) 

ofstream  readme ("taperchain. readme") ;  //  global  scope  fcts 

ofstream  Energylmpulse ( "taperchain. Enelmp") ; 

void  radii  (double  rlocalG)  { 
rlocal [nptles-1]  =  rlarge; 

if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles-1;  i++) 
rlocal [i]  =  rlarge; 

else 

for  (int  i  =  2;  i  <  nptles+1;  i++) 

rlocal [nptles-i]  =  rlocal [nptles-i+1]  *  (1  -  q*0.01); 
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> 


void  masses  (double  r  []  ,  double  masslocal  [])  { 
const  double  pi  =  4  *  atan(l.O); 

const  double  masslarge  =  (4. 0/3.0)  *  pi  *  rlarge*rlarge*rlarge  *  rho ; 
masslocal [nptles-1]  =  masslarge; 

if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles-1;  i++) 
masslocal [i]  =  masslarge; 

else 

for  (int  i  =  0;  i  <  nptles-1;  i++) 

masslocal [i]  =  r [i] *r [i] *r  [i]  *  masslarge  /  (rlarge*rlarge*rlarge) ; 

> 

void  strenghtfac  (double  r  []  ,  double  alocal  [])  { 
alocal[0]  =  (2.0  /  (5.0  *  D))  *  (sqrt (r  [0] ) )  ; 
alocal [nptles]  =  (2.0  /  (5.0  *  D))  *  (sqrt (r [nptles-1] ))  ; 
if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  1;  i  <  nptles;  i++) 

alocal [i]  =  (2.0  /  (5.0*D))  *  (sqrt (0 . 5*rlarge) ) ; 

else 

for  (int  i  =  1;  i  <  nptles;  i++) 

alocal [i]  =  (2.0  /  (5.0  *  D))  *  (sqrt ( (r [i] *r [i-1] )/(r  [i] +r  [i-1] ) ) ) ; 

> 

//  initialpos  prints  absolute  initial  positions,  not  for  calculations 
void  initialpos  (double  r [] ,  double  xlnitiallocal [] )  { 

if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  (2.0*(i+l)  -  1)  *  rlarge; 
else  { 

xlnitiallocal [0]  =  r [0]  ; 

for  (int  i  =  1;  i  <  nptles;  i++) 
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xlnitiallocal [i]  =  xlnitiallocal [i-1]  +  r[i-l]  +  r[i]; 

> 

> 

//  absolutpos  prints  absolute  positions  to  ptle  files,  not  for  calculations 
void  absolutpos  (double  r [] ,  double  x [] , 
double  xlnitial[],  double  xAbsolutlocal  [] )  { 
for  (int  i  =  0;  i  <  nptles;  i++) 

xAbsolutlocal [i]  =  xlnitial[i]  +  x [i] ; 

} 

void  computeAccelerations  (double  x [] ,  double  a[],  double  r [] , 
double  acc[],  double  overbef ore [] , 
double  mass[] ,  double&  pot)  { 

pot  =0.0;  //  every  call  calculates  new  pot  contributions 

for  (int  i  =  0;  i  <  nptles;  i++)  //  zeroing  all  acc  in  every  call 

acc[i]  =0.0;  //  (=  every  timestep) 

/*******  potential/force  between  neighboring  ptles  *****************/ 
for  (int  i  =  0;  i  <  nptles-1;  i++)  { 

if  (x[i]  >  x[i+l])  {  //  only  when  overlap 

double  over  =  x[i]  -  x[i+l]; 
double  overnml  =  pow(over,  (xn  -  1.0)); 
pot  +=  over  *  overnml  *  a[i+l]; 
double  forceBetw  =  a[i+l]  *  xn  *  overnml; 

double  forceFactor; 

if  (overbef ore  [i+1]  <  over)  //  when  compressing 

forceFactor  =  1.0; 

else  forceFactor  =  epsilon;  //  when  decompressing 
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forceBetw  *=  forceFactor; 


acc[i]  -=  forceBetw; 
acc[i+l]  +=  forceBetw; 

overbef ore [i+1]  =  over; 

} 

else  overbef ore  [i+1]  =  0.0; 


//  dim  acc:  force 

//  sign(-) :  towards  smaller  x 

//  sign(+) :  towards  larger  x 

//  update  for  next  timestep 

//  reset  when  no  overlap 


/**  potential/force  between  fixed  wall  (small,  x=0)  <->  small  ptle  **/ 
if  (x  [0]  <  0)  { 

double  over  =  -  x [0] ; 

double  overnml  =  pow(over,  (xn  -  1.0)); 

pot  +=  over  *  overnml  *  a[0]; 

double  forceSmall  =  a[0]  *  xn  *  overnml; 


double  forceFactor; 
if  (overbef ore  [0]  <  over) 
forceFactor  =  1.0; 
else  forceFactor  =  epsilon; 

forceSmall  *=  forceFactor; 

acc[0]  +=  forceSmall; 

overbef ore  [0]  =  over; 

> 

else  overbef ore  [0]  =  0.0; 

/***  potential/force  between  fixed  wall  (large)  <->  large  ptle  * *****/ 
if  (x[nptles-l]  >  0)  { 

double  over  =  x[nptles-l]; 
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double  overnml  =  pow(over,  (xn  -  1.0)); 

pot  +=  over  *  overnml  *  a[nptles]; 

double  forceLarge  =  a[nptles]  *  xn  *  overnml; 

double  forceFactor; 
if  (overbef ore [nptles]  <  over) 
forceFactor  =  1.0; 
else  forceFactor  =  epsilon; 

forceLarge  *=  forceFactor; 

acc [nptles-1]  -=  forceLarge; 

overbef ore [nptles]  =  over; 

> 

else  overbef ore [nptles]  =  0.0; 


/*****  real  dim  of  acc:  division  by  mass  **********/ 
for  (int  i  =  0;  i  <  nptles;  i++) 
acc[i]  /=  mass[i]; 

} 

void  velocityVerletStep  (double  x [] ,  double  v[],  double  acc [] , 
double  a[],  double  r  []  ,  double  overbef  ore  []  , 
double  mass  [] ,  double&  pot)  { 

for  (int  j  =  0;  j  <  nptles;  j++)  { 

x[j]  +=  v[j]  *  dt  +  0.5  *  acc[j]  *  dt*dt; 
v[j]  +=  0.5  *  acc[j]  *  dt ; 

> 

computeAccelerations  (x,  a,  r,  acc,  overbefore,  mass,  pot); 
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for  (int  j  =  0;  j  <  nptles;  j++) 
v[j]  +=  0.5  *  acc[j]  *  dt ; 

> 

void  ptleHeader  (of stream*  print,  int  k)  { 

(*  print)  <<  "#  ptle  "  <<  k+1  <<  time  (musec)"  <<  ’ \t ’  <<  "x  (mm)" 

<<  ’ \t ’  <<  "v  (mm/musec) "  <<  1 \t 1  <<  "a  (mm/musec~2) " 

«  ’ \t ’  «  "kin.  E.  (J)"  «  ’ \t ’  «  "f  (kN)"  «  ’ \t ’ 

<<  "impulse  (mg*mm/musec) "  <<  1 \t ’  «  "xRelative  (mm)"  <<  ’\n’; 

} 

void  dumpData  (double  t,  double  mass[],  double  v [] ,  double  acc[], 
double  r [] ,  double  x [] ,  //  scope  absolutpos 

double  xlnitial [] ,  double  xAbsolut [] , 
of stream*  print,  int  k)  { 

double  keDumpPt [nptles] ;  //  new  arrays  for  dumping  data 

double  vDumpPt [nptles] ;  //  since  arrays  pass  by  argument 

double  accDumpPt [nptles] ;  //  dump  data  manipulated 

double  xDumpPt [nptles]  ; 

keDumpPt  [k]  =0.5*  mass  [k]  *  v [k] *v  [k] ; 

if  (keDumpPt [k]  <  1.0e-20)  keDumpPt [k]  =0.0;  //  set  small  values  to  zero 

vDumpPt [k]  =  v [k] ; 

if  (vDumpPt [k]  <  1.0e-20  kk  vDumpPt [k]  >  -1.0e-20)  vDumpPt [k]  =  0.0; 

accDumpPt  [k]  =  acc  [k] ; 

if  (accDumpPt [k]  <  1.0e-20  kk  accDumpPt  [k]  >  -1.0e-20)  accDumpPt [k]  =  0.0 

xDumpPt  [k]  =  x [k] ; 

if  (xDumpPt [k]  <  1.0e-20  kk  xDumpPt [k]  >  -1.0e-20)  xDumpPt [k]  =  0.0; 
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absolutpos  (r,  x,  xlnitial,  xAbsolut) ;  //  calculate  absolute  pos. 


(*  print)  <<  t  <<  ’\t’  <<  xAbsolut [k]  <<  ’ \t ’  «  vDumpPt  [k]  <<  ’ \t; 

<<  accDumpPt  [k]  <<  ’ \t ’  <<  keDumpPt  [k]  <<  ’ \t; 

<<  accDumpPt [k] *mass [k]  <<  ’ \t; 

<<  mass [k] * vDumpPt  [k]  <<  1 \t 1  <<  xDumpPt [k]  <<  ’  \n’; 

> 

void  dumpEnergylmpulse  (double  t,  double  kelocal,  double  telocal,  double  pot, 
double  ptotallocal,  double  mass  [] ,  double  v  [] )  { 
kelocal  =  0.0; 
ptotallocal  =  0.0; 

double  absptotallocal  =0.0;  //  scope  only  within  function 

for  (int  j  =  0;  j  <  nptles;  j++)  { 
kelocal  +=  mass[j]  *  v[j]*v[j]; 
ptotallocal  +=  mass[j]  *  v [ j ] ; 
absptotallocal  +=  mass[j]  *  abs(v[j]); 

> 

kelocal  *=  0.5; 
telocal  =  kelocal  +  pot; 
double  potDump  =  pot ; 

if  (kelocal  <  1.0e-20)  kelocal  =0.0;  //  set  very  small  values  to  zero 

if  (ptotallocal  <  1.0e-20  &&  ptotallocal  >  -1.0e-20)  ptotallocal  =  0.0; 
if  (telocal  <  1.0e-20)  telocal  =  0.0; 
if  (potDump  <  1.0e-20)  potDump  =  0.0; 
if  (absptotallocal  <  1.0e-20)  absptotallocal  =  0.0; 

Energy Impulse .precision(16) ; 

Energy Impulse  <<  t  <<  1 \t ’  «  kelocal  <<  1 \t ’  «  potDump  <<  1 \t ’ 

«  telocal  <<  ’ \t ’  <<  absptotallocal  <<  ’ \t; 

<<  ptotallocal  <<  ’\n’; 

> 

88 


void  readmelnfo  (double  ke ,  double  pot,  double  te,  double  ptotal, 
double  kelin,  double  kenin,  double  plin,  double  pnin, 
double  r  []  ,  double  xlnitial[],  double  mass  []  , 
double  a[])  { 


readme 

« 

’  \t 1  «  ***  TAPERCHAIN 

within  walls 

*** 

(-:" 

«  ’\n’ 

« 

’  \n’  ; 

readme 

« 

"parameter  of  this  run:  "  << 

’  \n’ 

« 

’  \n’ 

i 

readme 

« 

"total  number  of  particles: 

"  « 

’Xt 

’  « 

’Xt’ 

<<  nptles  << 

’Xn 

readme 

« 

"density  of  particles  (mg/mm 

'3)  : 

"  «  >\t 

’  « 

rho  <<  ’ \n’ ; 

readme 

« 

"quantity  D  of  particles  (mm 

~2/N) 

.  II 

«  ’ 

Xt’ 

«  D  «  ’  \n’  ; 

readme 

« 

"radius  of  large  ptle  (mm) : 

"  « 

’Xt 

’  « 

’Xt’ 

<<  rlarge  << 

’Xn 

readme 

« 

"tapering  factor  (°/0)  :  "  <<  1 

\t’  «  ’ 

\t’  «  ’\t’  «  q  «  ’ \n 

J  . 
i 

readme 

« 

"exponent  in  potential:  "  << 

’Xt’ 

« 

’Xt’ 

« 

’\t’  «  xn  « 

’Xn 

readme 

« 

"timestepwidth  (musec) :  "  << 

’\t’ 

« 

’Xt’ 

« 

’\t’  «  dt  « 

’Xn 

readme 

« 

"#  steps  integration  loop:  " 

«  1 

Xt’ 

«  ’ 

Xt’ 

<<  nsteps  <<  ’ 

Xn’ 

readme 

« 

"stepwidth  diagnostics:  "  << 

’Xt' 

« 

’Xt’ 

« 

’ \t ’  <<  "every  " 

<<  idiagp  <<  "  timesteps"  <<  ’\n’; 


readme 

« 

"stepwidth  dump:  "  << 

C+ 

A 

A 

/ 

c+ 

« 

’Xt’ 

<<  "every  " 

<<  idump  << 

"  timesteps"  <<  ’\n’ 

f 

readme 

« 

"initial  v  small  ptle 

(mm/musec) : 

"  « 

’Xt’ 

«  vlin  « 

’Xn’  ; 

readme 

« 

"initial  v  large  ptle 

(mm/musec) : 

"  « 

’Xt’ 

<<  vnin  << 

’Xn’  ; 

readme 

« 

"restitution  factor  for  all  ptles: 

"  « 

’Xt 

’  <<  1-epsilon  <<  ’ \n 

«  ’\n’; 

readme 

« 

"total  length  of  run 

(musec) :  "  << 

’Xt’ 

« 

dt  *  nsteps 

«  ’\n’; 

readme 

« 

"total  rows  recorded 

for  .Enelmp  file: 

"  « 

’\t’  «  ’\t 

J 

<<  nsteps/idiagp  +1  <<  ’\n’; 

readme 

« 

"total  rows  recorded 

for  particle 

files 

.  II 

«  ’\t’ 

<<  nsteps/idump  +1  <<  ’ \n ’  «  ’\n’; 
readme  <<  "Initial  system  info  (t=0) :  "  <<  ’\n’  <<  ’\n’; 
readme  <<  "kin.  E.  (J)"  <<  ’\t’  <<  "pot.  E.  (J)"  <<  1 \t ’ 

<<  "tot.  E.  (J)"  <<  ’\t’  <<  "total  impulse  (mg*mm/musec) "  <<  ’\n’; 
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readme  <<  ke  <<  ’ \t ’  <<  pot  <<  1 \t ’  <<  te  <<  1 \t ’  <<  ptotal 
«  »\n’  «  ’ \n’  ; 

readme  <<  "kin.  E.  of  small  particle  (J) :  "  <<  ’\t’  <<  ’ \t ’  <<  kelin 
«  ’  \n’  ; 

readme  <<  "kin.  E.  of  large  particle  (J) :  "  <<  ’\t’  <<  ’ \t ’  <<  kenin 
«  ’  \n’  ; 

readme  <<  "impulse  of  small  particle  (mg*mm/musec) :  "  <<  ’\t’ 

<<  plin  <<  ’ \n’  ; 

readme  <<  "impulse  of  large  particle  (mg*mm/musec) :  "  <<  ’\t’ 

<<  pnin  <<  ’ \n ’  <<  ’\n’; 
readme  <<  "particle  radii  (mm):  "  <<  ’\n’; 
for  (int  i=0;  i  <  nptles;  i++)  { 
readme  <<  r[i]  <<  ’ \tJ; 

> 

readme  <<  ’ \n ’  <<  ’ \n’ ; 

readme  <<  "initial  particle  positions  (mm):  "  <<  ’\n’; 
for  (int  i=0;  i  <  nptles;  i++)  { 
readme  <<  xlnitial[i]  <<  1 \t ’ ; 

> 

readme  <<  ’ \n ’  <<  ’ \n’ ; 

readme  <<  "total  length  of  one  dimensional  alignment  (mm):  "  <<  ’\t; 

<<  xlnitial [nptles-1]  +  r[nptles-l]  <<  ’ \n ’  <<  ’\n’; 
readme  <<  "particle  masses  (mg):  "  <<  ’Xn’; 
for  (int  i=0;  i  <  nptles;  i++)  { 
readme  <<  mass[i]  <<  1 \t ’  ; 

> 

readme  <<  ’ \n ’  <<  ’ \n’ ; 

readme  <<  "particle  interaction  strenghts  (0 . 0316*N/mm~ (3/2) ) :  "  <<  ’\nJ 
for  (int  i=0;  i  <  nptles+1;  i++)  { 
readme  <<  a[i]  <<  ’Xt1; 

> 

readme  <<  ’ \n’; 


> 
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int  main  (  )  { 


double  r[nptles],  x[nptles],  xAbsolut [nptles]  ,  xlnitial [nptles]  , 
v[nptles],  acc [nptles],  mass [nptles]  ,  a[nptles+l], 
overbef ore [nptles+1] ,  pot; 

/********************  functions  *****************************/ 

radii  (r) ; 

masses  (r,  mass); 

strenghtfac  (r,  a); 

initialpos  (r,  xlnitial); 

/**********************  output  files  ****************************/ 

/******  ptle-f iles  *********/ 

ofstrecLm  print  [nptles]  ; 

for  (int  k  =  0;  k  <  nptles;  k++)  { 

string  filename; 
ostringstream  buffer; 

buffer  <<  "taperchain_"  <<  k+1  <<  ".dat"; 
filename  =  buffer . str() ; 

print [k] . open(f ilename . c_str() ) ;  //  convert  string  to  char 

> 

for  (int  k  =  0;  k  <  nptles;  k++)  //  header  for  particles 

ptleHeader(&  print [k] ,  k) ; 
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//  header  for  . Enelmp  file 

Energy Impulse  <<  "#  time"  <<  ’ \t;  <<  "kin.  E.  (J)"  <<  1 \t ’ 

«  >\t;  «  "pot.  E.  (J)"  «  ’ \t ’  «  ’ \t ’ 

«  "total  E.  (J)"  «  ’\t;  «  >\t> 

«  " I (total  imp.) I  (mg*mm/musec) "  <<  ’ \tJ 
<<  "total  imp.  (mg*mm/musec) "  <<  ’\n’; 

for  (int  i  =  0;  i  <  nptles;  i++)  {  //  zeroing 

x[i]  =0.0;  //  relative  particle  positions  for  calculation 

v[i]  =  0.0; 
acc[i]  =  0.0; 
overbef ore  [i]  =  0.0; 

> 

overbef ore [nptles]  =  0.0; 
double  t  =  0.0; 

v[0]  =  vlin;  //  mind  special  input  data 

v[nptles-l]  =  vnin; 

/****************  initial  energy  info  *********************************/ 
double  kelin  =  0.5  *  mass[0]  *  v [0] *v  [0] ;  //  initial  kE  of  edge  ptles 

double  kenin  =  0.5  *  mass [nptles-1]  *  v [nptles-1]  *v [nptles-1]  ; 
double  plin  =  mass[0]  *  v  [0] ;  //  initial  impulse  of  edge  particles 

double  pnin  =  mass [nptles-1]  *  v [nptles-1]; 

computeAccelerations  (x,  a,  r,  acc,  overbefore,  mass,  pot); 

//  call  for  initial  potential  energy 

double  ke  =  0.0;  //  checking  for  initial  system  energy 

double  ptotal  =  0.0; 

for  (int  i  =  0;  i  <  nptles;  i++)  { 

acc[i]  =0.0;  //  reset  acc  for  verlet  in  case  of  initial  overlap 
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ke  +=  mass[i]  *  v[i]*v[i]; 
ptotal  +=  mass[i]  *  v [i] ; 

> 

ke  *=  0.5; 

double  te  =  pot  +  ke ;  //  initial  total  system  energy 

readmelnfo  (ke,  pot,  te,  ptotal,  kelin,  kenin,  plin,  pnin, 
r,  xlnitial,  mass,  a); 

/*******************  begin  of  timestep  loop  *********************/ 
for  (unsigned  int  i  =  0;  i  <  nsteps+1;  i++)  { 

t  =  i  *  dt ; 

velocityVerletStep  (x,  v,  acc,  a,  r,  overbefore,  mass,  pot); 

/***********  check  system  energy,  plot  data  *******************/ 
if  ((i  7,  idiagp)  ==  0)  { 

dumpEnergyImpulse(t ,  ke,  te,  pot,  ptotal,  mass,  v) ; 

} 

/*********  check  particle  energy,  plot  data  *******************/ 
if  ((i  7.  idump)  ==  0)  { 

for  (int  k  =  0;  k  <  nptles;  k++)  { 

dumpData(t,  mass,  v,  acc,  //  scope  dumpData 

r,  x,  xlnitial,  xAbsolut ,  //  scope  absolutpos 

&  print [k] ,  k) ; 

> 

} 

}-/********************  end  of  timestep  loop  *********************/ 

readme . close () ; 

Energy Impulse . close () ; 
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for  (int  k  =  0;  k  <  nptles;  k++)  //  closing  particle  files 


print [k] . close ()  ; 


> 
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Appendix  B 


DTC  Code  Modifications 


/**************************  ALTERABLE  PARAMETER:  ********* *******************/ 
int  nptles  =3;  //  total  number  of  particles 

const  double  f  =  0.9;  //  fractional  size  of  interstitial  grain  w.r.t  last  grain 

const  double  rho=4.42;  //  TiAlV  (mg/mm~3) 
const  double  D  =  0.01206;  //  TiAlV  (mnT2/N) 

const  double  rlarge  =5.0;  //  (radius  of  large  ptle  (mm)) 

const  double  q  =  0.0;  //  (tapering  factor  (%)) 

const  double  xn  =  2.5;  //  (exponent  in  potential) 

const  double  dt  =  0.00001;  //  (timestepwidth  (musec)) 

const  unsigned  int  nsteps  =  100000000;  //  (#  steps  integration  loop) 
const  int  idiagp  =  20000;  //  (stepwidth  diagnostics) 

const  int  idump  =  20000;  //  (stepwidth.  dump) 

const  double  vlin  =0.0;  //  (initial  v  small  ptle  (mm/musec)) 

const  double  vnin  =  -0.01;  //  initial  v  large  ptle  (mm/musec)) 

const  double  epsilon  =1.0;  //  ((1  -  restitution  factor)  all  ptles) 

ofstream  readme ("taperchain. readme") ;  //  global  scope  fcts 

ofstream  Energylmpulse ( "taperchain. Enelmp") ; 

//  Generate  radii  and  masses  for  DTC 

void  spheres  (double  rlocal[],  double  masslocalG)  { 

rlocal [nptles-1]  =  rlarge;  //  shifts  everything  to  index  starting  at  zero 
double  tapering  =  1  -  q*0.01; 
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const  double  pi  =  4  *  atan(l.O); 

const  double  masslarge  =  (4. 0/3.0)  *  pi  *  pow(rlarge,3)  *  rho; 
masslocal [nptles-1]  =  masslarge; 
if  (q==0  kk  f==1.0)  //  Monodisperse  in  DTC 
for  (int  i=0;  i<nptles-l;  i++)  -[ 
rlocal[i]  =  rlarge; 
masslocal [i]  =  masslarge; 

} 

else  if  (q==0)  {  //  Quasi-Monodisperse :  Avoid  roundoff  errors  without  tapering 
for  (int  i=l;  i<nptles-l;  i=i+2)  {  //  Interstitial  grains 
rlocal[i]  =  rlarge*f ; 

masslocal[i]  =  (4. 0/3.0)  *  pi  *  powCrlocal [i] ,3)  *  rho; 

> 

for  (int  i=0;  i<nptles-l;  i=i+2)  {  //  Non-interstitial  grains 
rlocal[i]  =  rlarge; 
masslocal [i]  =  masslarge; 

> 

> 

//  non-monodisperse  chains 
else  { 

for  (int  i=l;  i<nptles-l;  i=i+2)  {  //  Find  radii  of  interstitial  grains 
rlocal[i]  =  f*pow (tapering, (nptles-l)/2)*rlarge; 
masslocal[i]  =  (4. 0/3.0)  *  pi  *  pow(rlocal [i] ,3)  *  rho; 

> 

for  (int  i=nptles-l;  i>=0;  i=i-2)  {  //Find  radii  of  non-interstitital  grains 
rlocal[i-2]  =  rlocal[i]  *  tapering; 
masslocal [i— 2]  =  (4. 0/3.0)  *  pi  *  pow(rlocal [i-2] ,3)  *  rho; 

> 

> 

> 

void  strenghtfac  (double  r [] ,  double  alocal[])  { 
alocal[0]  =  (2.0  /  (5.0  *  D))  *  (sqrt (r  [0] ) )  ; 
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alocal [nptles]  =  (2.0  /  (5.0  *  D))  *  (sqrt (r [nptles-1] ) ) ; 
if  (q  ==  0  &&  f  ==  1 . 0)  //  avoid  roundoff  errors  w/out  tapering 
for  (int  i  =  1;  i  <  nptles;  i++) 

alocal [i]  =  (2.0  /  (5.0*D))  *  (sqrt (0 . 5*rlarge) ) ; 

else 

for  (int  i  =  1;  i  <  nptles;  i++) 

alocal [i]  =  (2.0  /  (5.0  *  D))  *  (sqrt ( (r [i] *r [i-1] )/(r  [i] +r  [i-1] ) ) ) ; 

} 

//  initialpos  prints  absolute  initial  positions,  not  for  calculations 
void  initialpos  (double  r [] ,  double  xlnitiallocal [] )  { 

if  (q  ==  0  &&  f  ==  1 . 0  )  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  (2.0*(i+l)  -  1)  *  rlarge; 
else  { 

xlnitiallocal [0]  =  r [0] ; 

for  (int  i  =  1;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  xlnitiallocal [i-1]  +  r[i-l]  +  r  [i] ; 

> 

> 
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Appendix  C 

PERL  Code 


# ! /usr/bin/env  perl 

#  This  program  automates  the  considerable  task  of  setting  up  parametric  studies  on  the  tapered 

#  spherical  elastic  1-D  grain  problem.  It  takes  an  input  file:  taperchain29 . cpp  and  searches 

#  through  it  replacing  the  values  of  epsilon,  q,  and  N  in  for  loops  and  spitting  out  a  file 

#  in  the  appropriate  directory.  The  directories  are  created  on  the  fly. 

#  This  version  uses  the  updated  directories  and  is  looking  at  initial  velocity  on  the  small 

#  grain. 

$w  =  0.1; 

$FILE_NAME  =  "taperchain28_w01 . cpp" ; 

$S0URCE_DIR  =  "/home/rldoney" ; 

#$S0URCE_DIR  =  "/Users/bob/Work/Classes/Spring  2004/Dr.  Sen  study/runs"; 

$FILE_IN  =  "$SOURCE_DIR/$FILE_NAME" ; 

$0UT_DIR  =  "/nf s/scratch/rldoney/TiAlV/D . SimpleTapered" ; 

#$0UT_DIR  =  "$SOURCE_DIR/TiAlV/D . SimpleTapered/D . Vin_large/" ; 

$EPSIL0N_CHK  =  "const  double  epsilon  =";  #  Set  pattern  to  match  line  with  epsilon 
$Q_CHK  =  "const  double  q  =" ;  #  Set  pattern  to  match  line  with  q 

$N_CHK  =  "int  nptles=";  #  Set  pattern  to  match  line  with  N 

systemC'clear") ;  #  clear  the  screen; 

#for($w  =  0.0;  $w<=0.02;  $w+=0.01)  {  #  Restitution 

$epsilon  =  1.0  -  $w;  #  taperchain . cpp  program  uses  epsilon  instead  of  w  directly 
print "\n\n: :::::::  w=$w  \t  epsilon  =  $epsilon  :::::::: :\n"; 
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system ( " date  ;  +DATE :  Zm/Zd/ZyZnTIME :  %H :  7.M :  %S  ’  "  )  ; 
print "\n" ; 

chdir ( " $OUT_DIR" )  or  die  "Cant  open  $0UT_DIR" ; 
mkdir("w$w") ; 

chdir("w$w")  or  die  "Cant  open  $0UT_DIR/w$w" ; 

for($N=3;  $N  <=20;  $N++)  { 
mkdir("N$N") ; 

chdir("N$N")  or  die  "Cant  open  N$N" ; 

for($q=0;  $q  <=10;  $q++)  { 
mkdir("N$N\q$q") ; 

chdir("N$N\q$q")  or  die  "Cant  open  N$N\q$q" ; 

$CURRENT_DIR  =  ‘pwd‘ ; 
chop($CURRENT_DIR) ;  #  remove  trailing  \n 
print"Current  Dir:  $CURRENT_DIR\n" ; 

$FILE_0UT=  "$CURRENT_DIR/$FILE_NAME" ; 

open (FROM,  "$FILE_IN")  or  die  "Cant  open  $FILE_IN:  $!"; 

open(T0 ,  ">$FILE_0UT" )  or  die  "Cant  open  $CURRENT_DIR/$FILE_OUT:  $!"; 

print"  "; 

system  ("date  ;+TIME:  °/„H :  °/„M :  °/„S  ’  "  )  ; 

while (<FR0M>)  {  #  Read  in  the  file  line  by  line  into  $_ 

#  Replacements  (Regular  expression  matching) 
if($w  !=  0)  { 

s/$EPSIL0N_CHK  \d+\ . \d+/$EPSIL0N_CHK  $epsilon/;  #  change  epsilon 

> 

s/$N_CHK\d+/$N_CHK  $N/;  #  change  N 
s/$Q_CHK  \d+\ . \d+/$Q_CHK  $q\.0/;  #  change  q 

print  TO  #  Write  the  current  line  with  any  changes  to  TO 
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close  TO; 


} 

close  FROM; 

#  Escape  to  the  shell,  compile  the  file,  and  run  it 

#  It  was  unexpected,  but  we  need  the  ’  ’  because  of  the  spaces  in  the  path 
system  ("g++  1 $FILE_0UT’ ") ; 

system  ("./a. out"); 
chdir(" 

} 

chdir (" 
print ("\n") ; 

> 

system ( " date  ;  +DATE :  7.m/7.d/7.y%nTIME :  7.H :  %M :  7.S  ’  "  )  ; 
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Appendix  D 

Energy  partitioning  in  the  DTC 
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t=10us 


Figure  D.l:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  10/rs  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.2:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  20fj,s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.3:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  30/j.s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.4:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  40/j.s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.5:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  50/j.s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.6:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  60/j.s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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t=70us 


Figure  D.7:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  70fxs  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 


108 


t=80[.is 


1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 

grain  number 


1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 


Figure  D.8:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  80fj,s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.9:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  90fi,s  where 
qd  =  {0,0.05,0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Figure  D.10:  Instantaneous  kinetic  energy  per  grain  for  various  DTC  configurations  at  t  =  100/xs 
where  qd  =  {0,  0.05,  0.1}  and  /  =  {1.0,  0.7,  0.3}. 
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Appendix  E 

Normalization 


Normalization  has  posed  some  challenges  in  trying  to  properly  assess  the  absorption  quality  of  a 
simple  tapered  chain  (STC).  Proper  normalization  schemes  will  help  the  architect  better  determine 
which  STC  is  best  for  which  application.  Recall  that  the  goal  is  to  measure  the  energy  at  the  last 
grain  versus  the  energy  put  into  the  system  by  the  first.  In  general,  the  functional  form  will  stay  the 
same  regardless  of  the  strategy  and  adjustments  in  the  normalization  will  simply  scale  the  kinetic 
energy  (KE)  surface.  In  some  cases,  the  output  force  of  each  chain  is  based  on  the  maximum  value 
felt  by  any  possible  chain  under  consideration  (i.e.  monodisperse  and  no  energy  loss)1.  This  gives  a 
measure  of  how  one  chain  is  better  than  another. 

In  this  communication,  we  have  chosen  to  form  the  ratio  based  on  the  output  KE  and  force 
felt  by  each  specific  chain.  This  serves  to  grade  the  individual  effectiveness  of  any  chain  without 
reference  to  another.  In  choosing  a  peak  value  with  this  method,  one  could  identify  either  when 
the  impulse  first  hits  the  last  grain  or  look  for  the  absolute  maximum  peak  whenever  it  may  occur. 
We  have  chosen  to  use  the  former  for  several  reasons.  First,  it  ignores  the  complexity  of  nonlinear 
reverberations  which  can  lead  to  large  peaks  at  unpredictable  times.  Second,  we  argue  that  this  is 
just  as  realistic  as  selecting  the  maximum  value  anywhere  in  the  time  spectrum. 

It  turns  out  that  in  most  cases,  the  absolute  maximum  is  the  first  peak.  There  are  special  cases 
where  the  maximum  may  occur  at  later  times  and  this  needs  to  be  investigated  further.  In  figure  2.5 
(in  the  body  of  the  report)  for  q  =  0.1,  for  example,  we  see  the  striking  occurrence  of  the  secondary 
pulse  about  225  fj,s  being  much  stronger  than  the  initial  arrival  at  35  fis.  This  is  one  of  those 
instances  that  disagrees  with  the  way  we  choose  to  normalize  our  KE  surfaces.  We  have  investigated 
this  particular  case  further  without  including  extraneous  plots  and  report  the  following  observations. 
The  effect  exists  for  N  =  15  —  20  for  constant  to.  When  N  =  20  is  held  fixed  and  restitution  is 
increased,  the  peak  KE  once  again  occurs  for  the  first  arrival  of  the  pulse.  As  q  increases,  so  do  the 
number  of  collisions  and  the  requisite  energy  loss  (since  u>  ^  0).  Therefore,  it  is  less  likely  to  find  a 

global  peak  later  in  the  simulation.  The  situation  is  further  complicated  by  the  interplay  between  q 
1  Pfannes,  J.  Energy  propagation  in  granular  chains  M.S.  Thesis,  SUNY  Buffalo,  May  2003. 
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Appendix  F 


CUBIT  script  for  3D  STC 


#  This  is  a  cubit  script  to  generate  a  Simple 

#  Tapered  Chain  (STC) . 

#  units  are  in  mm 

#{scalef  =  1.0e-3} 

#{N  =  5}  ##  Number  of  grains 

#{n  =  N-l}  ##  used  only  for  loop  control 

#{nn  =  N+l}  ##  For  vertices 

#{Ri  =5.0}  ##  Largest  grain  radius 

#{q  =  0.07}  ##  tapering  (percentage) 

#{nr  =  6}  ##  No.  of  radial  divisions 

#{ntheta  =  40}  ##  No.  of  azimuthal  divisions  (must  be  multiple  of  8) 

#{wall_length  =  Ri*2}  ##  Length  of  wall  boundary  along  x 

#{wall_y  =  Ri*2}  ##  Transverse  y  size  of  wall 

#{wall_z  =  wall_y}  ##  Transverse  z  size  of  wall 

graphics  mode  wire 
create  sphere  radius  {Ri} 

#{rold  =  Ri}  ##  Set  the  initial  r_old 

#{xold  =  0}  ##  Initial  center  of  1st  grain 

#{i  =  1}  ##  counter  for  cubit  body  number 

-(loop(n)}  ##  can’t  loop  over  algebraic  ops 

#{i  =  i  +1} 

#{rnew  =  rold  *  (1.0-q)}  ##  determine  size  of  next  grain 
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#{xnew  =  xold  +  rold  +  mew} 
create  sphere  radius  {rnew} 
body  {i}  move  {xnew}  0  0 
#{rold  =  rnew} 

#{xold  =  xnew} 

{EndLoop} 

#{edge=  xnew  +  rnew} 

#create  vertex  {edge}  0  0 

create  brick  x  {wall_length}  y 
body  {nn}  move  {wall_length/2 . 
body  {nn}  size  1.0 
mesh  volume  {nn} 

#merge  all 


##  find  position  of  next  grain 

##  move  the  grain  along  x 
##  update  values  and  prepare 
##  for  next  new  values 

##  position  of  far  side  of  chain 

{wall_y}  z  {wall_z} 

)  +  edge}  0  0 


#{i  =  0} 

{loop(N)} 

#{i  =i+l} 

volume  {i}  scheme  sphere  graded_interval  {nr}  Az_interval  {ntheta}  fraction  0.7 
mesh  volume  {i} 

volume  {i}  smooth  scheme  laplacian 
smooth  volume  {i} 

{EndLoop} 

##  Scale  - 

Transform  Mesh  Output  scale  {scalef} 

##  BLOCK  ASSIGNMENT:  spheres  +  block  - 

#{i  =  0} 

{loop(nn)} 

#{i  =i+l} 

block  {i}  volume  {i} 
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{EndLoop} 


###  BCs  - 

#{i  =  0} 

{loop(N)} 

#{i  =i+l} 

nodeset  {i}0  surface  {i} 
{EndLoop} 


#{i  =  0> 

{loop(N)} 

#{i  =i+l} 

sideset  {i}0  surface  {i} 

{EndLoop} 

##  Block-sphere  interface:  that  surface  always  seems  to  be  N+4 
nodeset  999  surface  {N+4} 
sideset  999  surface  {N+4} 

mesh  visibility  off 
disp 

export  genesis  ’stc.gen’  overwrite 
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Appendix  G 


ALEGRA  script  for  3D  STC 


$$#  Alegra  input  script  for  STC  dynamics 

$$# 

$$#  Bob  Doney 
$$#  10.29.2004 
$$# 

$$#  05.30.2006  updated  to  3d  lagrange 

$$# 

$$# - 


$$#  -  Grain  setup - 

${scalef  =  1.0e-3}  ##  Scale  SI  to  mm 

${N  =  5}  ##  Number  of  grains 

${n  =  N-l> 

${nn  =  N+l} 

${Ri  =  5.0*scalef}  ##  Largest  grain  radius 

${q  =  0.07}  ##  tapering  (percentage) 


title:  STC 
units,  si 

termination  cycle,  1 
$$#termination  time,  1000. 0e-6 

solid  dynamics 

cartesian  3d 
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ignore  kinematic  errors 


domain 

remesh  iterations  10 
remap  iterations  1 
mhis  advection 
smyra  interface  tracker 
$$#  slic  interface  tracker 

initial  remesh  movement  limiter  0.005 
remesh  movement  ratio  1.025 
remesh  movement  limiter  0.4 
end 

$$#  -  Attempt  to  use  boundaries  as  walls  that  constrain  the  chain 

no  displacement,  sideset  999,  x 

$$#  GLOBAL  CONTACT  -  using  SNL  input  to  start  with 

global  contact  $$#  for  3D  applications 

{loop(N)} 

${i=j*10} 
sideset  {i} 

j=j+l} 

-CEndLoop} 

sideset  999 

package  =  acme 

search  algorithm  =  augmented  dynamic  search 
enforcement  algorithm  =  td  enforcement 
enforcement  iterations  =  20 
default  data 

kinematic  partition  =  auto 
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friction  model  id  =  1 


search  normal  tolerance  =  .0001 
search  tangential  tolerance  =  .0001 
end 

friction  model  1  frictionless 
$no  subkeywords 
end 
end 

initial  block  velocity:  block  1,  x  500.0 
maximum  initial  time  step  1.0e-9 


$$#  TRACER  PARTICLES 
tracer  points 

lagrangian  tracer  1  x=0  y=0  z=0  $$#  Initial  one  at  origin 

${rold  =  Ri}  $$#  Set  the  initial  r_old 

${xold  =  0}  $$#  Initial  center  of  1st  grain 

${i=0} 

$$#  Center  of  spheres 
{loop(n)}  $$#  N-l 

${i  =  i+1} 

${rnew  =  rold  *  (1.0-q)}  $$#  determine  radius  of  next  grain 

${xnew  =  xold  +  rold  +  rnew} 

lagrangian  tracer  {i+1}  x  =  {xnew}  y=0  z=0 

${rold  =  rnew}  $$#  update  values  and  prepare 

${xold  =  xnew}  $$#  for  next  new  values 

{EndLoop} 

$$#  Wall-last  sphere  interface 
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lagrangian  tracer  {nil}  x=  {xnew+rnew}  y=0  z=0 


end 


$$#  End  tracer  section 


$$#  MESH 
${i=0> 

{loop(nn)} 

${i=i+l} 
block  {i} 

LAGRANGIAN  MESH 
material  {i} 

#  hourglass  control  {i} 
end 

{EndLoop} 


end  $$#  —  end  physics  — 

$$# - 

$$#  —  PLOTTING  — 

$$#double  precision  exodus 
emit  output,  time  interval=l . e-6 
emit  plot,  time  interval=l . e-6 
emit  hisplt,  time  interval=0 . 2e-6 

plot  variables 
no  underscores 
velocity 
density 
density, avg 
temperature 
temperature,  avg 
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energy 


stress 
stress,  avg 
end 

history  plot  variables 
no  material  globals 
end 

$$# - 

$$#  —  MATERIALS 


${i=0> 

{loop(nn)} 

${i=i+l} 

material  {i}  $$#  Copper  spheres 

model  =  100 
end 

{EndLoop} 

$—  ELASTIC-PLASTIC  — 
model  100  cth  elastic  plastic 
eos  model  =  11 

yield  model  =  12 

end 

$$#  EOS 

model  11,  keos  sesame 
feos  =  ’sesame’ 
neos  =  3320 
clip  =  1.0 
end 
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$$#  YIELD 

model  12  johnson  cook  ep 
matlabel  =  ’COPPER’ 
end 

exit 
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